++Figure 3. The graph shows the number of accumulated secondary +extinctions that occur when removing species from the most to the least +connected species +
+ +| 1 | -- |
- #' Extinctions analysis for ecological networks- |
-
| 2 | -- |
- #'- |
-
| 3 | -- |
- #' The SimulateExtinctions function, can be used to test how the order of species- |
-
| 4 | -- |
- #' extinctions, species-dependency on existing interaction strength, and rewiring potential might affect the stability of the network by comparing The extinction history- |
-
| 5 | -- |
- #' and checking for secondary extinctions.- |
-
| 6 | -- |
- #'- |
-
| 7 | -- |
- #' @param Network a network representation as a an adjacency matrix, edgelist,- |
-
| 8 | -- |
- #' or a network object- |
-
| 9 | -- |
- #' @param Method a character with the options Mostconnected and Ordered- |
-
| 10 | -- |
- #' @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.- |
-
| 11 | -- |
- #' @param NetworkType a character with the options Trophic and Mutualistic - is used to calculate secondary extinctions.- |
-
| 12 | -- |
- #' @param clust.method a character with the options cluster_edge_betweenness,- |
-
| 13 | -- |
- #' cluster_label_prop or cluster_infomap, defaults to cluster_infomap- |
-
| 14 | -- |
- #' @param IS either numeric or a named vector of numerics. Identifies the threshold of relative interaction strength which species require to not be considered secondarily extinct (i.e. IS = 0.3 leads to removal of all nodes which lose 70 percent of their interaction strength in the Network argument). If a named vector, names must correspond to vertex names in Network argument.- |
-
| 15 | -- |
- #' @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.- |
-
| 16 | -- |
- #' @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. If Rewiring == function(x){x}, this matrix is expected to contain probabilities of a connection being present between species-pairs.
- |
-
| 17 | -- |
- #' @param RewiringProb a numeric which identifies the threshold at which to assume rewiring potential is met.- |
-
| 18 | -- |
- #' @param verbose Logical. Whether to report on function progress or not.- |
-
| 19 | -- |
- #' @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.- |
-
| 20 | -- |
- #' @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. This process is repeated until the entire network has gone extinct.- |
-
| 21 | -- |
- #'- |
-
| 22 | -- |
- #' 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.- |
-
| 23 | -- |
- #'- |
-
| 24 | -- |
- #' 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.- |
-
| 25 | -- |
- #'- |
-
| 26 | -- |
- #' When clust.method = cluster_edge_betweenness computes the network modularity using cluster_edge_betweenness methods from igraph to detect communities- |
-
| 27 | -- |
- #' When clust.method = cluster_label_prop computes the network modularity using cluster_label_prop methods from igraph to detect communities- |
-
| 28 | -- |
- #' 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- |
-
| 29 | -- |
- #' @examples- |
-
| 30 | -- |
- #' # Mostconnected example- |
-
| 31 | -- |
- #' data("net")
- |
-
| 32 | -- |
- #' SimulateExtinctions(Network = net, Method = "Mostconnected",- |
-
| 33 | -- |
- #' clust.method = "cluster_infomap")- |
-
| 34 | -- |
- #'- |
-
| 35 | -- |
- #' #first Ordered example- |
-
| 36 | -- |
- #' data("net")
- |
-
| 37 | -- |
- #' SimulateExtinctions(Network = net, Order = c(1,2,3,4,5,6,7,8,9,10),- |
-
| 38 | -- |
- #' Method = "Ordered" , clust.method = "cluster_infomap")- |
-
| 39 | -- |
- #'- |
-
| 40 | -- |
- #' #Second Ordered example- |
-
| 41 | -- |
- #' data("net")
- |
-
| 42 | -- |
- #' SimulateExtinctions(Network = net, Order = c(2,8,9),- |
-
| 43 | -- |
- #' Method = "Ordered", clust.method = "cluster_infomap")- |
-
| 44 | -- |
- #'- |
-
| 45 | -- |
- #' #Network-Dependency Example- |
-
| 46 | -- |
- #' data("net")
- |
-
| 47 | -- |
- #' SimulateExtinctions(Network = net, Order = c(2,8), IS = 0.3,- |
-
| 48 | -- |
- #' Method = "Ordered", clust.method = "cluster_infomap")- |
-
| 49 | -- |
- #'- |
-
| 50 | -- |
- #' #Rewiring- |
-
| 51 | -- |
- #' data("net")
- |
-
| 52 | -- |
- #' data(dist)- |
-
| 53 | -- |
- #' SimulateExtinctions(Network = net, Order = c(2,8), IS = 0.3,- |
-
| 54 | -- |
- #' # assuming an exponential decline in rewiring potential- |
-
| 55 | -- |
- #' # as values in RewiringDist increase- |
-
| 56 | -- |
- #' Rewiring = function(x){1-pexp(x, rate = 1/0.5)},
- |
-
| 57 | -- |
- #' RewiringDist = dist, # distance matrix- |
-
| 58 | -- |
- #' RewiringProb = 0.2, # low threshold for rewiring potential- |
-
| 59 | -- |
- #' Method = "Ordered", clust.method = "cluster_infomap")- |
-
| 60 | -- |
- #'- |
-
| 61 | -- |
- #' #Rewiring, assuming dist contains probabilities- |
-
| 62 | -- |
- #' #' data("net")
- |
-
| 63 | -- |
- #' data(dist)- |
-
| 64 | -- |
- #' SimulateExtinctions(Network = net, Order = c(2,8), IS = 0.3,- |
-
| 65 | -- |
- #' Rewiring = function(x){x}, # no changes to the RewiringDist object means
- |
-
| 66 | -- |
- #' RewiringDist = dist, RewiringProb = 0.2,- |
-
| 67 | -- |
- #' Method = "Ordered", clust.method = "cluster_infomap")- |
-
| 68 | -- |
- #' @importFrom dplyr desc- |
-
| 69 | -- |
- #' @author Derek Corcoran <derek.corcoran.barrios@gmail.com>- |
-
| 70 | -- |
- #' @author M. Isidora Ávila-Thieme <msavila@uc.cl>- |
-
| 71 | -- |
- #' @author Erik Kusch <erik.kusch@bio.au.dk>- |
-
| 72 | -- |
- #' @export- |
-
| 73 | -- |
- SimulateExtinctions <- function(Network, Method, Order = NULL,- |
-
| 74 | -- |
- NetworkType = "Trophic", clust.method = "cluster_infomap",- |
-
| 75 | -- |
- IS = 0,- |
-
| 76 | -- |
- Rewiring = FALSE, RewiringDist, RewiringProb = 0.5,- |
-
| 77 | -- |
- verbose = TRUE){
- |
-
| 78 | -8x | -
- Network <- .DataInit(x = Network)- |
-
| 79 | -! | -
- if(!NetworkType %in% c("Trophic", "Mutualistic")){stop("Please specify NetworkType as either 'Trophic' or 'Mutualistic'")}
- |
-
| 80 | -- | - - | -
| 81 | -3x | -
- if(!is.null(Order)){Method <- "Ordered"}
- |
-
| 82 | -- | - - | -
| 83 | -8x | -
- '%ni%'<- Negate('%in%')
- |
-
| 84 | -! | -
- if(Method %ni% c("Mostconnected", "Ordered")) stop('Choose the right method. See ?SimulateExtinction.')
- |
-
| 85 | -- | - - | -
| 86 | -8x | -
- edgelist <- network::as.matrix.network.edgelist(Network,matrix.type="edgelist") #Prey - Predator- |
-
| 87 | -8x | -
- if(Method == "Mostconnected"){
- |
-
| 88 | -- |
- # if(NetworkType == "Trophic"){
- |
-
| 89 | -- |
- # Conected <- as.numeric(names(sort(table(edgelist[,1]), decreasing = TRUE)))- |
-
| 90 | -- |
- # }else{
- |
-
| 91 | -5x | -
- Grado <- NULL- |
-
| 92 | -5x | -
- Conected <- data.frame(ID = 1:network::network.size(Network), Grado = sna::degree(edgelist, c("total")))
- |
-
| 93 | -5x | -
- Conected <- dplyr::arrange(Conected, desc(Grado))$ID- |
-
| 94 | -- |
- # }- |
-
| 95 | -5x | -
- DF <- ExtinctionOrder(Network = Network, Order = Conected, clust.method = clust.method,- |
-
| 96 | -5x | -
- IS = IS, Rewiring = Rewiring, RewiringDist = RewiringDist,- |
-
| 97 | -5x | -
- verbose = verbose, RewiringProb = RewiringProb, NetworkType = NetworkType,- |
-
| 98 | -5x | -
- RecalcConnect = TRUE)- |
-
| 99 | -- |
- }- |
-
| 100 | -8x | -
- if(Method == "Ordered"){
- |
-
| 101 | -3x | -
- DF <- ExtinctionOrder(Network = Network, Order = Order, clust.method = clust.method,- |
-
| 102 | -3x | -
- IS = IS, Rewiring = Rewiring, RewiringDist = RewiringDist,- |
-
| 103 | -3x | -
- verbose = verbose, RewiringProb = RewiringProb, NetworkType = NetworkType)- |
-
| 104 | -- |
- }- |
-
| 105 | -- | - - | -
| 106 | -8x | -
- return(DF)- |
-
| 107 | -- |
- }- |
-
| 108 | -- | - - | -
| 109 | -- |
- #' Extinctions analysis from custom order- |
-
| 110 | -- |
- #'- |
-
| 111 | -- |
- #' 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).- |
-
| 112 | -- |
- #'- |
-
| 113 | -- |
- #' @param Network a network representation as a an adjacency matrix, edgelist, or a network object- |
-
| 114 | -- |
- #' @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.- |
-
| 115 | -- |
- #' @param NetworkType a character with the options Trophic and Mutualistic - is used to calculate secondary extinctions.- |
-
| 116 | -- |
- #' @param clust.method a character with the options cluster_edge_betweenness,- |
-
| 117 | -- |
- #' cluster_label_prop or cluster_infomap, defaults to cluster_infomap- |
-
| 118 | -- |
- #' @param IS either numeric or a named vector of numerics. Identifies the threshold of relative interaction strength which species require to not be considered secondarily extinct (i.e. IS = 0.3 leads to removal of all nodes which lose 70percent of their interaction strength in the Network argument). If a named vector, names must correspond to vertex names in Network argument.- |
-
| 119 | -- |
- #' @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.- |
-
| 120 | -- |
- #' @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. If Rewiring == function(x){x}, this matrix is expected to contain probabilities of a connection being present between species-pairs.
- |
-
| 121 | -- |
- #' @param RewiringProb a numeric which identifies the threshold at which to assume rewiring potential is met.- |
-
| 122 | -- |
- #' @param verbose Logical. Whether to report on function progress or not.- |
-
| 123 | -- |
- #' @param RecalcConnect Logical. Whether to recalculate connectedness of each node following each round of extinction simulation and subsequently update extinction order with newly mostconnected nodes.- |
-
| 124 | -- |
- #' @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.- |
-
| 125 | -- |
- #' @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.- |
-
| 126 | -- |
- #'- |
-
| 127 | -- |
- #' When clust.method = cluster_edge_betweenness computes the network modularity using cluster_edge_betweenness methods from igraph to detect communities- |
-
| 128 | -- |
- #' When clust.method = cluster_label_prop computes the network modularity using cluster_label_prop methods from igraph to detect communities- |
-
| 129 | -- |
- #' 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- |
-
| 130 | -- |
- #'- |
-
| 131 | -- |
- #' @importFrom dplyr relocate- |
-
| 132 | -- |
- #' @importFrom ggplot2 aes_string- |
-
| 133 | -- |
- #' @importFrom ggplot2 geom_line- |
-
| 134 | -- |
- #' @importFrom ggplot2 ggplot- |
-
| 135 | -- |
- #' @importFrom ggplot2 theme_bw- |
-
| 136 | -- |
- #' @importFrom ggplot2 xlab- |
-
| 137 | -- |
- #' @importFrom ggplot2 ylab- |
-
| 138 | -- |
- #' @importFrom network as.matrix.network.edgelist- |
-
| 139 | -- |
- #' @importFrom network delete.vertices- |
-
| 140 | -- |
- #' @importFrom network network.edgecount- |
-
| 141 | -- |
- #' @importFrom network network.size- |
-
| 142 | -- |
- #' @importFrom network network.density- |
-
| 143 | -- |
- #' @importFrom network get.vertex.attribute- |
-
| 144 | -- |
- #' @importFrom network get.edge.attribute- |
-
| 145 | -- |
- #' @importFrom igraph as.undirected- |
-
| 146 | -- |
- #' @importFrom igraph E- |
-
| 147 | -- |
- #' @importFrom sna degree- |
-
| 148 | -- |
- #' @importFrom stats complete.cases- |
-
| 149 | -- |
- #' @importFrom network as.matrix.network.adjacency- |
-
| 150 | -- |
- #' @importFrom igraph graph_from_adjacency_matrix- |
-
| 151 | -- |
- #' @importFrom igraph cluster_edge_betweenness- |
-
| 152 | -- |
- #' @importFrom igraph cluster_label_prop- |
-
| 153 | -- |
- #' @importFrom igraph cluster_infomap- |
-
| 154 | -- |
- #' @importFrom igraph modularity- |
-
| 155 | -- |
- #' @importFrom stats na.omit- |
-
| 156 | -- |
- #' @importFrom utils setTxtProgressBar- |
-
| 157 | -- |
- #' @importFrom utils txtProgressBar- |
-
| 158 | -- |
- #' @importFrom dplyr desc- |
-
| 159 | -- |
- #' @author Derek Corcoran <derek.corcoran.barrios@gmail.com>- |
-
| 160 | -- |
- #' @author M. Isidora Ávila-Thieme <msavila@uc.cl>- |
-
| 161 | -- |
- #' @author Erik Kusch <erik.kusch@bio.au.dk>- |
-
| 162 | -- |
- #' @export- |
-
| 163 | -- |
- ExtinctionOrder <- function(Network, Order, NetworkType = "Trophic", clust.method = "cluster_infomap",- |
-
| 164 | -- |
- IS = 0,- |
-
| 165 | -- |
- Rewiring = FALSE, RewiringDist, RewiringProb = 0.5,- |
-
| 166 | -- |
- verbose = TRUE,- |
-
| 167 | -- |
- RecalcConnect = FALSE- |
-
| 168 | -- |
- ){
- |
-
| 169 | -! | -
- if(!NetworkType %in% c("Trophic", "Mutualistic")){stop("Please specify NetworkType as either 'Trophic' or 'Mutualistic'")}
- |
-
| 170 | -- |
- # Setting up Objects for function run ++++++++++ ++++++++++ ++++++++++ ++++++++++- |
-
| 171 | -129x | -
- Link_density <- Modularity <- Grado <- NULL- |
-
| 172 | -129x | -
- Network <- .DataInit(x = Network)- |
-
| 173 | -129x | -
- edgelist <- network::as.matrix.network.edgelist(Network,matrix.type="edgelist") #Prey - Predator- |
-
| 174 | -129x | -
- Conected <- data.frame(ID = 1:network::network.size(Network), Grado = sna::degree(edgelist, c("total")))
- |
-
| 175 | -129x | -
- Conected <- dplyr::arrange(Conected, desc(Grado))- |
-
| 176 | -129x | -
- Conected1 <- Order- |
-
| 177 | -- | - - | -
| 178 | -- |
- ## Interaction Strength Loss Preparations ++++++++++ ++++++++++- |
-
| 179 | -129x | -
- if(length(IS )== 1){ # when the same dependency is to be applied for all species
- |
-
| 180 | -129x | -
- IS <- rep(IS, network::network.size(Network)) # repeat IS argument for all species- |
-
| 181 | -129x | -
- names(IS) <- network::get.vertex.attribute(Network, "vertex.names") # assign species names to IS arguments- |
-
| 182 | -- |
- }else{
- |
-
| 183 | -! | -
- 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.")}
- |
-
| 184 | -- |
- }- |
-
| 185 | -- |
- ## Rewiring Preparations ++++++++++ ++++++++++- |
-
| 186 | -129x | -
- if(!isFALSE(Rewiring)){ # if rewiring has been specified
- |
-
| 187 | -! | -
- if(!exists("RewiringDist")){stop("To execute rewiring simulations, you need to specify the RewiringDist argument as well as the Rewiring argument.")}
- |
-
| 188 | -1x | -
- diag(RewiringDist)<- NA # set diagonal to NA as one can't rewire to oneself- |
-
| 189 | -1x | -
- if(length(Rewiring) == 1){ # when the same Rewiring happen for all species
- |
-
| 190 | -1x | -
- fun <- deparse1(Rewiring) # turn function into string- |
-
| 191 | -1x | -
- Rewiring <- rep(fun, network::network.size(Network)) # repeat function for each species- |
-
| 192 | -1x | -
- names(Rewiring) <- network::get.vertex.attribute(Network, "vertex.names") # assign names of species in network to rewiring functions- |
-
| 193 | -- |
- }- |
-
| 194 | -! | -
- 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.")}
- |
-
| 195 | -- |
- }- |
-
| 196 | -- | - - | -
| 197 | -- |
- # Base net calculations ++++++++++ ++++++++++ ++++++++++ ++++++++++- |
-
| 198 | -- |
- ## base interaction strengths per node ++++++++++ ++++++++++- |
-
| 199 | -129x | -
- options(warn=-1) # turn warning off- |
-
| 200 | -129x | -
- Weight_mat <- net <- network::as.matrix.network.adjacency(Network, attrname = "weight")- |
-
| 201 | -129x | -
- options(warn=0) # turn warnings on- |
-
| 202 | -129x | -
- if(sum(IS) != 0){
- |
-
| 203 | -- |
- # if(sum(network::get.edge.attribute(Network, "weight"), na.rm = TRUE) == 0){
- |
-
| 204 | -- |
- # 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.")
- |
-
| 205 | -- |
- # }- |
-
| 206 | -1x | -
- netgraph <- suppressMessages(igraph::graph_from_adjacency_matrix(net, weighted = TRUE))- |
-
| 207 | -1x | -
- if(NetworkType == "Trophic"){
- |
-
| 208 | -1x | -
- strengthbaseout <- igraph::strength(netgraph, mode = "out")- |
-
| 209 | -1x | -
- strengthbasein <- igraph::strength(netgraph, mode = "in")- |
-
| 210 | -- |
- }- |
-
| 211 | -1x | -
- if(NetworkType == "Mutualistic"){
- |
-
| 212 | -! | -
- strengthbasenet <- igraph::strength(netgraph)- |
-
| 213 | -- |
- }- |
-
| 214 | -- |
- }- |
-
| 215 | -- | - - | -
| 216 | -- |
- ## identification of producers and top predators ++++++++++ ++++++++++- |
-
| 217 | -129x | -
- indegreebasenet <- sna::degree(Network, cmode = "indegree")- |
-
| 218 | -129x | -
- indegreebasenetzeros <- sum(sna::degree(Network, cmode = "indegree") == 0)- |
-
| 219 | -129x | -
- indegreetopnetzeros <- sum(sna::degree(Network, cmode = "outdegree") == 0)- |
-
| 220 | -129x | -
- Producers <- network::get.vertex.attribute(Network, "vertex.names")[sna::degree(Network, cmode = "indegree") == 0]- |
-
| 221 | -129x | -
- TopPredators <- network::get.vertex.attribute(Network, "vertex.names")[sna::degree(Network, cmode = "outdegree") == 0]- |
-
| 222 | -- | - - | -
| 223 | -- |
- ## output object ++++++++++ ++++++++++- |
-
| 224 | -129x | -
- DF <- data.frame(Spp = rep(NA, length(Order)),- |
-
| 225 | -129x | -
- S = rep(NA, length(Order)),- |
-
| 226 | -129x | -
- L = rep(NA, length(Order)),- |
-
| 227 | -129x | -
- C = rep(NA, length(Order)),- |
-
| 228 | -129x | -
- Link_density = rep(NA, length(Order)),- |
-
| 229 | -129x | -
- SecExt = rep(NA,length(Order)),- |
-
| 230 | -129x | -
- Pred_release = rep(NA,length(Order)),- |
-
| 231 | -129x | -
- Iso_nodes = rep (NA,length(Order)),- |
-
| 232 | -129x | -
- Modularity = rep (NA,length(Order)))- |
-
| 233 | -129x | -
- Secundaryext <- c()- |
-
| 234 | -129x | -
- Predationrel <- c()- |
-
| 235 | -129x | -
- accExt <- c()- |
-
| 236 | -129x | -
- totalExt <- c()- |
-
| 237 | -129x | -
- FinalExt <- list()- |
-
| 238 | -129x | -
- Conected3 <- c()- |
-
| 239 | -- | - - | -
| 240 | -- |
- # Sequential extinction simulation ++++++++++ ++++++++++ ++++++++++ ++++++++++- |
-
| 241 | -9x | -
- if(verbose){ProgBar <- txtProgressBar(max = length(Order), style = 3)}
- |
-
| 242 | -129x | -
- for (i in 1:length(Order)){
- |
-
| 243 | -- |
- # print(i)- |
-
| 244 | -- |
- ### creating temporary network + deleting vertices if they have been set to go extinct ++++++++++ ++++++++++- |
-
| 245 | -3299x | -
- if(length(accExt)==0){ # on first iteration
- |
-
| 246 | -129x | -
- Temp <- Network- |
-
| 247 | -129x | -
- DF$Spp[i] <- Conected1[i]- |
-
| 248 | -129x | -
- network::delete.vertices(Temp, c(DF$Spp[1:i]))- |
-
| 249 | -- |
- }- |
-
| 250 | -3299x | -
- if (length(accExt)>0){ # on any subsequent iteration
- |
-
| 251 | -3170x | -
- Temp <- Network- |
-
| 252 | -3170x | -
- Temp <- network::delete.vertices(Temp, c(accExt))- |
-
| 253 | -3170x | -
- edgelist <- network::as.matrix.network.edgelist(Temp,matrix.type="edgelist")- |
-
| 254 | -- | - - | -
| 255 | -3170x | -
- if(RecalcConnect){
- |
-
| 256 | -87x | -
- Conected2 <- data.frame(ID = 1:network::network.size(Temp), Grado = sna::degree(edgelist, c("total")))
- |
-
| 257 | -87x | -
- Conected2 <- arrange(Conected2, desc(Grado))- |
-
| 258 | -87x | -
- for(j in sort(accExt)){
- |
-
| 259 | -1673x | -
- Conected2$ID <- ifelse(Conected2$ID < j, Conected2$ID, Conected2$ID + 1)- |
-
| 260 | -- |
- }- |
-
| 261 | -87x | -
- DF$Spp[i] <- Conected2$ID[1]- |
-
| 262 | -- |
- }else{
- |
-
| 263 | -3083x | -
- DF$Spp[i] <- Conected1[i]- |
-
| 264 | -- |
- }- |
-
| 265 | -- | - - | -
| 266 | -3170x | -
- Temp <- Network- |
-
| 267 | -3170x | -
- network::delete.vertices(Temp, unique(c(c(DF$Spp[1:i]),accExt)))- |
-
| 268 | -- |
- }- |
-
| 269 | -- | - - | -
| 270 | -- |
- ## network metrics to output object ++++++++++ ++++++++++- |
-
| 271 | -3299x | -
- DF$S[i] <- network::network.size(Temp)- |
-
| 272 | -3299x | -
- DF$L[i] <- network::network.edgecount(Temp)- |
-
| 273 | -3299x | -
- DF$C[i] <- network::network.density(Temp)- |
-
| 274 | -3299x | -
- DF$Link_density[i] <- DF$L[i]/DF$S[i]- |
-
| 275 | -- | - - | -
| 276 | -- |
- ## premature complete annihilation message ++++++++++ ++++++++++- |
-
| 277 | -3299x | -
- if(i > 1){
- |
-
| 278 | -3170x | -
- if(DF$L[i-1] == 0){
- |
-
| 279 | -9x | -
- if(verbose){setTxtProgressBar(ProgBar, length(Order))}
- |
-
| 280 | -127x | -
- warning(paste("Your network become completely unconnected before all primary extinctions were simulated. This happened at extinction step", i-1, "out of", length(Order)))
- |
-
| 281 | -127x | -
- break- |
-
| 282 | -- |
- }- |
-
| 283 | -- |
- }- |
-
| 284 | -- | - - | -
| 285 | -- |
- ## calculating modularity ++++++++++ ++++++++++- |
-
| 286 | -3172x | -
- Networkclass = class(Temp)- |
-
| 287 | -3172x | -
- if (Networkclass[1] == "matrix"){
- |
-
| 288 | -! | -
- netgraph = igraph::graph_from_adjacency_matrix(Temp, mode = "directed", weighted = TRUE)- |
-
| 289 | -- |
- }- |
-
| 290 | -3172x | -
- if (Networkclass[1] == "network"){
- |
-
| 291 | -3172x | -
- net = network::as.matrix.network.adjacency(Temp)- |
-
| 292 | -3172x | -
- netgraph = suppressMessages(igraph::graph_from_adjacency_matrix(net, mode = "directed", weighted = TRUE))- |
-
| 293 | -- |
- }- |
-
| 294 | -3172x | -
- if (clust.method == "cluster_edge_betweenness"){
- |
-
| 295 | -5x | -
- Membership = suppressWarnings(cluster_edge_betweenness(netgraph, weights = igraph::E(- |
-
| 296 | -5x | -
- igraph::as.undirected(netgraph)- |
-
| 297 | -5x | -
- )$weight, directed = TRUE, edge.betweenness = TRUE,- |
-
| 298 | -5x | -
- merges = TRUE, bridges = TRUE, modularity = TRUE, membership = TRUE))- |
-
| 299 | -3167x | -
- }else if (clust.method == "cluster_label_prop"){
- |
-
| 300 | -5x | -
- Membership = suppressWarnings(cluster_label_prop(netgraph, weights = igraph::E(igraph::as.undirected(netgraph))$weight, initial = NULL,- |
-
| 301 | -5x | -
- fixed = NULL))- |
-
| 302 | -3162x | -
- }else if (clust.method == "cluster_infomap"){
- |
-
| 303 | -3162x | -
- nb.trials = network::network.size(Temp)- |
-
| 304 | -3162x | -
- Membership = suppressWarnings(igraph::cluster_infomap(igraph::as.undirected(netgraph),- |
-
| 305 | -3162x | -
- e.weights = igraph::E(- |
-
| 306 | -3162x | -
- igraph::as.undirected(netgraph)- |
-
| 307 | -3162x | -
- )$weight,- |
-
| 308 | -3162x | -
- v.weights = NULL,- |
-
| 309 | -3162x | -
- nb.trials = nb.trials,- |
-
| 310 | -3162x | -
- modularity = TRUE))- |
-
| 311 | -- | - - | -
| 312 | -! | -
- } else if (clust.method == "none"){
- |
-
| 313 | -! | -
- Membership = NA- |
-
| 314 | -! | -
- }else stop('Select a valid method for clustering. ?SimulateExtinction')
- |
-
| 315 | -3172x | -
- DF$Modularity[i] <- Membership$modularity- |
-
| 316 | -- | - - | -
| 317 | -- |
- ## rewiring ++++++++++ ++++++++++- |
-
| 318 | -3172x | -
- accExt <- unique(append(accExt, DF$Spp[1:i]))- |
-
| 319 | -3172x | -
- if(!isFALSE(Rewiring)){
- |
-
| 320 | -- |
- ### identify rewiring potential ++++++++++- |
-
| 321 | -36x | -
- Rewiring_df <- data.frame(Direction = NA,- |
-
| 322 | -36x | -
- Species = NA,- |
-
| 323 | -36x | -
- NewPartner = NA,- |
-
| 324 | -36x | -
- LostPartner = NA,- |
-
| 325 | -36x | -
- IS = NA)- |
-
| 326 | -36x | -
- Rewiring_df <- na.omit(Rewiring_df)- |
-
| 327 | -- |
- #### loop over all deleted vertices and the connections lost because of their exclusion- |
-
| 328 | -36x | -
- for(Iter_PrimaryExt in 1:length(accExt)){
- |
-
| 329 | -- |
- # Iter_PrimaryExt = 1- |
-
| 330 | -794x | -
- LostPartner <- network::get.vertex.attribute(Network, "vertex.names")[accExt[Iter_PrimaryExt]] # name of primary extinction species- |
-
| 331 | -794x | -
- LostISCol <- Weight_mat[, LostPartner] # lost interaction strength with nodes now slated for secondary extinction- |
-
| 332 | -794x | -
- LostISRow <- Weight_mat[LostPartner, ]- |
-
| 333 | -794x | -
- Lost_df <- data.frame(LostIS = c(LostISCol, LostISRow),- |
-
| 334 | -794x | -
- Direction = rep(c(1,2), c(length(LostISCol), length(LostISRow))),- |
-
| 335 | -794x | -
- names = c(names(LostISCol), names(LostISRow))- |
-
| 336 | -- |
- )- |
-
| 337 | -794x | -
- Lost_df <- Lost_df[Lost_df$LostIS != 0, ]- |
-
| 338 | -794x | -
- if(nrow(Lost_df)!=0){
- |
-
| 339 | -794x | -
- for(Iter_LostIS in 1:nrow(Lost_df)){ ## looping over all species that were linked to the current primary extinction
- |
-
| 340 | -- |
- # Iter_LostIS = 1- |
-
| 341 | -22661x | -
- if(Rewiring[which(names(Rewiring) == Lost_df$names[Iter_LostIS])] == "function (x) { x }"){
- |
-
| 342 | -22661x | -
- LostPartnerSim <- eval(str2lang(Rewiring[which(names(Rewiring) == Lost_df$names[Iter_LostIS])]))(RewiringDist[,Lost_df$names[Iter_LostIS]]) # probability of rewiring to each node- |
-
| 343 | -- |
- }else{
- |
-
| 344 | -! | -
- LostPartnerSim <- eval(str2lang(Rewiring[which(names(Rewiring) == Lost_df$names[Iter_LostIS])]))(RewiringDist[,LostPartner]) # probability of rewiring to each node in network given rewiring function and species similraity- |
-
| 345 | -- |
- }- |
-
| 346 | -22661x | -
- names(LostPartnerSim) <- colnames(RewiringDist)- |
-
| 347 | -22661x | -
- 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- |
-
| 348 | -22661x | -
- RewiredPartner <- names(which.max(RewiringCandidates)) # most likely rewiring partner- |
-
| 349 | -22661x | -
- if(!is.null(RewiredPartner)){ # if a rewired partner has been found
- |
-
| 350 | -523x | -
- Rewiring_df <- rbind(Rewiring_df,- |
-
| 351 | -523x | -
- data.frame(Direction = Lost_df$Direction[Iter_LostIS],- |
-
| 352 | -523x | -
- Species = Lost_df$names[Iter_LostIS],- |
-
| 353 | -523x | -
- NewPartner = RewiredPartner,- |
-
| 354 | -523x | -
- LostPartner = LostPartner,- |
-
| 355 | -523x | -
- IS = Lost_df$LostIS[Iter_LostIS])- |
-
| 356 | -- |
- )- |
-
| 357 | -- |
- }- |
-
| 358 | -- |
- }- |
-
| 359 | -- |
- }- |
-
| 360 | -- |
- }- |
-
| 361 | -- | - - | -
| 362 | -- |
- ### shift rewired interaction strengths ++++++++++- |
-
| 363 | -36x | -
- if(nrow(Rewiring_df) != 0){
- |
-
| 364 | -- |
- #### shift interaction weights in Weight_mat- |
-
| 365 | -34x | -
- for(Iter_Rewiring in 1:nrow(Rewiring_df)){
- |
-
| 366 | -- |
- # Iter_Rewiring = 1- |
-
| 367 | -- |
- ## assigning shifted interaction strength- |
-
| 368 | -523x | -
- ColSpec <- Rewiring_df[Iter_Rewiring,4-Rewiring_df[Iter_Rewiring,"Direction"]]- |
-
| 369 | -523x | -
- RowSpec <- Rewiring_df[Iter_Rewiring,1+Rewiring_df[Iter_Rewiring,"Direction"]]- |
-
| 370 | -523x | -
- Weight_mat[RowSpec, ColSpec] <- Weight_mat[RowSpec, ColSpec] + Rewiring_df[Iter_Rewiring,"IS"]- |
-
| 371 | -- |
- ## deleting shiften interaction strength- |
-
| 372 | -- | - - | -
| 373 | -523x | -
- ColLost <- ifelse(Rewiring_df[Iter_Rewiring, "Direction"] == 1,- |
-
| 374 | -523x | -
- Rewiring_df[Iter_Rewiring, "LostPartner"],- |
-
| 375 | -523x | -
- Rewiring_df[Iter_Rewiring, "Species"])- |
-
| 376 | -523x | -
- RowLost <- ifelse(Rewiring_df[Iter_Rewiring, "Direction"] == 1,- |
-
| 377 | -523x | -
- Rewiring_df[Iter_Rewiring, "Species"],- |
-
| 378 | -523x | -
- Rewiring_df[Iter_Rewiring, "LostPartner"])- |
-
| 379 | -523x | -
- Weight_mat[RowLost, ColLost] <- 0- |
-
| 380 | -- |
- }- |
-
| 381 | -- |
- #### establishing rewired network and deleting primary extinction nodes- |
-
| 382 | -34x | -
- Network <- as.network(Weight_mat, matrix.type = "adjacency", ignore.eval=FALSE, names.eval='weight')- |
-
| 383 | -34x | -
- Temp <- Network- |
-
| 384 | -34x | -
- network::delete.vertices(Temp, unique(c(c(DF$Spp[1:i]),accExt)))- |
-
| 385 | -- |
- }- |
-
| 386 | -- |
- }- |
-
| 387 | -- | - - | -
| 388 | -- |
- ## identify secondary extinctions ++++++++++ ++++++++++- |
-
| 389 | -- |
- ### Relative Interaction Strength loss ++++++++++- |
-
| 390 | -3172x | -
- if(sum(IS) == 0){
- |
-
| 391 | -3141x | -
- SecundaryextNames <- network::get.vertex.attribute(Temp, "vertex.names")[which(sna::degree(Temp) == 0)]- |
-
| 392 | -3141x | -
- Secundaryext <- match(SecundaryextNames, network::get.vertex.attribute(Network, "vertex.names"))- |
-
| 393 | -- |
- }else{
- |
-
| 394 | -31x | -
- if(NetworkType == "Trophic"){
- |
-
| 395 | -31x | -
- AbsISin <- igraph::strength(suppressMessages(igraph::graph_from_adjacency_matrix(- |
-
| 396 | -31x | -
- network::as.matrix.network.adjacency(Temp, attrname = "weight"),- |
-
| 397 | -31x | -
- weighted = TRUE)- |
-
| 398 | -31x | -
- ), mode = "in")- |
-
| 399 | -31x | -
- RelISRemainIn <- AbsISin / strengthbasein[names(strengthbasein) %in% network::get.vertex.attribute(Temp, "vertex.names")]- |
-
| 400 | -31x | -
- SecundaryextNames <- names(which(AbsISin == 0 | RelISRemainIn < IS[match(names(RelISRemainIn), names(IS))]))- |
-
| 401 | -31x | -
- SecundaryextNames <- SecundaryextNames[!(SecundaryextNames %in% Producers)]- |
-
| 402 | -31x | -
- Secundaryext <- match(SecundaryextNames, network::get.vertex.attribute(Network, "vertex.names"))- |
-
| 403 | -- |
- }- |
-
| 404 | -- | - - | -
| 405 | -31x | -
- if(NetworkType == "Mutualistic"){
- |
-
| 406 | -! | -
- AbsIS <- igraph::strength(suppressMessages(igraph::graph_from_adjacency_matrix(- |
-
| 407 | -! | -
- network::as.matrix.network.adjacency(Temp, attrname = "weight"),- |
-
| 408 | -! | -
- weighted = TRUE)- |
-
| 409 | -- |
- ))- |
-
| 410 | -! | -
- RelISloss <- AbsIS / strengthbasenet[names(strengthbasenet) %in% network::get.vertex.attribute(Temp, "vertex.names")]- |
-
| 411 | -! | -
- SecundaryextNames <- names(which(AbsIS == 0 | RelISloss < IS[match(names(RelISloss), names(IS))]))- |
-
| 412 | -! | -
- Secundaryext <- match(SecundaryextNames, network::get.vertex.attribute(Network, "vertex.names"))- |
-
| 413 | -- |
- }- |
-
| 414 | -- |
- }- |
-
| 415 | -- | - - | -
| 416 | -- |
- ### for trophic networks ++++++++++- |
-
| 417 | -3172x | -
- if(NetworkType == "Trophic"){
- |
-
| 418 | -3172x | -
- MidPredExt <- network::get.vertex.attribute(Temp, "vertex.names")[sna::degree(Temp, cmode = "indegree") == 0]- |
-
| 419 | -3172x | -
- MidPredExt <- match(MidPredExt[!(MidPredExt %in% Producers)], network::get.vertex.attribute(Network, "vertex.names"))- |
-
| 420 | -3172x | -
- SecundaryextTrue <- unique(c(SecundaryextNames[!(SecundaryextNames %in% as.character(Producers))],- |
-
| 421 | -3172x | -
- MidPredExt))- |
-
| 422 | -3172x | -
- Secundaryext <- match(SecundaryextTrue, network::get.vertex.attribute(Network, "vertex.names"))- |
-
| 423 | -3172x | -
- DF$SecExt[i] <- length(Secundaryext)- |
-
| 424 | -3172x | -
- DF$Pred_release[i] <- length(SecundaryextNames[!(SecundaryextNames %in% as.character(TopPredators))])- |
-
| 425 | -- |
- }- |
-
| 426 | -- |
- ### for mutualistic networks ++++++++++- |
-
| 427 | -3172x | -
- if(NetworkType == "Mutualistic"){
- |
-
| 428 | -! | -
- DF$SecExt[i] <- length(Secundaryext)- |
-
| 429 | -- |
- }- |
-
| 430 | -3172x | -
- DF$Iso_nodes[i] <- sum(sna::degree(Temp) == 0)- |
-
| 431 | -- | - - | -
| 432 | -- |
- ## Return of objects ++++++++++ ++++++++++- |
-
| 433 | -3172x | -
- FinalExt[[i]] <- Secundaryext- |
-
| 434 | -3172x | -
- accExt <- append(accExt, DF$Spp[1:i])- |
-
| 435 | -3172x | -
- accExt <- unique(append(accExt,Secundaryext))- |
-
| 436 | -107x | -
- if(verbose){setTxtProgressBar(ProgBar, i)}
- |
-
| 437 | -- | - - | -
| 438 | -- |
- ## final Temp deletion ++++++++++ ++++++++++- |
-
| 439 | -3172x | -
- if(length(accExt)>0){
- |
-
| 440 | -3172x | -
- Temp <- Network- |
-
| 441 | -3172x | -
- Temp <- network::delete.vertices(Temp, c(accExt))- |
-
| 442 | -3172x | -
- edgelist <- network::as.matrix.network.edgelist(Temp,matrix.type="edgelist")- |
-
| 443 | -- |
- # DF$Spp[i] <- Conected1[i]- |
-
| 444 | -3172x | -
- Temp <- Network- |
-
| 445 | -3172x | -
- network::delete.vertices(Temp, unique(c(c(DF$Spp[1:i]),accExt)))- |
-
| 446 | -- |
- }- |
-
| 447 | -- |
- }- |
-
| 448 | -- | - - | -
| 449 | -- |
- # return of final data objects ++++++++++ ++++++++++- |
-
| 450 | -129x | -
- DF <- DF[!is.na(DF$Spp),]- |
-
| 451 | -129x | -
- DF$AccSecExt <- cumsum(DF$SecExt)- |
-
| 452 | -129x | -
- DF$NumExt <- 1:nrow(DF)- |
-
| 453 | -129x | -
- DF$TotalExt <- DF$AccSecExt + DF$NumExt- |
-
| 454 | -129x | -
- DF <- relocate(DF, Modularity, .after = Link_density)- |
-
| 455 | -129x | -
- class(DF) <- c("data.frame", "SimulateExt")
- |
-
| 456 | -- | - - | -
| 457 | -129x | -
- return(list(sims = DF,- |
-
| 458 | -129x | -
- Network = Temp))- |
-
| 459 | -- |
- }- |
-
| 460 | -- | - - | -
| 461 | -- |
- #' Random extinction- |
-
| 462 | -- |
- #'- |
-
| 463 | -- |
- #' 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.- |
-
| 464 | -- |
- #'- |
-
| 465 | -- |
- #' @param Network a network representation as a an adjacency matrix, edgelist,- |
-
| 466 | -- |
- #' or a network object- |
-
| 467 | -- |
- #' @param nsim numeric, number of simulations- |
-
| 468 | -- |
- #' @param Record logical, if TRUE, records every simulation and you can read the- |
-
| 469 | -- |
- #' raw results in the object FullSims- |
-
| 470 | -- |
- #' @param plot logical if TRUE, will add a graph to the results- |
-
| 471 | -- |
- #' @param SimNum numeric, how many nodes to register for primary extinction. By default sets all of them.- |
-
| 472 | -- |
- #' @param NetworkType a character with the options Trophic and Mutualistic - is used to calculate secondary extinctions.- |
-
| 473 | -- |
- #' @param clust.method a character with the options cluster_edge_betweenness,- |
-
| 474 | -- |
- #' cluster_label_prop or cluster_infomap, defaults to cluster_infomap- |
-
| 475 | -- |
- #' @param parallel if TRUE, it will use parallel procesing, if FALSE (default) it will run- |
-
| 476 | -- |
- #' sequentially- |
-
| 477 | -- |
- #' @param ncores numeric, number of cores to use if using parallel procesing- |
-
| 478 | -- |
- #' @param IS either numeric or a named vector of numerics. Identifies the threshold of relative interaction strength which species require to not be considered secondarily extinct (i.e. IS = 0.3 leads to removal of all nodes which lose 70 precent of their interaction strength in the Network argument). If a named vector, names must correspond to vertex names in Network argument.- |
-
| 479 | -- |
- #' @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.- |
-
| 480 | -- |
- #' @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. If Rewiring == function(x){x}, this matrix is expected to contain probabilities of a connection being present between species-pairs.
- |
-
| 481 | -- |
- #' @param RewiringProb a numeric which identifies the threshold at which to assume rewiring potential is met.- |
-
| 482 | -- |
- #' @param verbose Logical. Whether to report on function progress or not.- |
-
| 483 | -- |
- #' @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 95percent 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.- |
-
| 484 | -- |
- #' @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.- |
-
| 485 | -- |
- #'- |
-
| 486 | -- |
- #' When clust.method = cluster_edge_betweenness computes the network modularity using cluster_edge_betweenness methods from igraph to detect communities- |
-
| 487 | -- |
- #' When clust.method = cluster_label_prop computes the network modularity using cluster_label_prop methods from igraph to detect communities- |
-
| 488 | -- |
- #' 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- |
-
| 489 | -- |
- #'- |
-
| 490 | -- |
- #' @examples- |
-
| 491 | -- |
- #' #first example- |
-
| 492 | -- |
- #' \dontrun{
- |
-
| 493 | -- |
- #' data("More_Connected")
- |
-
| 494 | -- |
- #' RandomExtinctions(Network = More_Connected, nsim = 20)- |
-
| 495 | -- |
- #'- |
-
| 496 | -- |
- #' # Using parallel procesing- |
-
| 497 | -- |
- #' ## Detect your number of cores divide by 2- |
-
| 498 | -- |
- #'- |
-
| 499 | -- |
- #' cores <- ceiling(parallel::detectCores()/2)- |
-
| 500 | -- |
- #'- |
-
| 501 | -- |
- #' RandomExtinctions(Network = More_Connected, nsim = 20, parallel = TRUE, ncores = cores)- |
-
| 502 | -- |
- #' }- |
-
| 503 | -- |
- #'- |
-
| 504 | -- |
- #' @importFrom doParallel registerDoParallel- |
-
| 505 | -- |
- #' @importFrom dplyr group_by- |
-
| 506 | -- |
- #' @importFrom dplyr mutate- |
-
| 507 | -- |
- #' @importFrom dplyr summarise- |
-
| 508 | -- |
- #' @importFrom foreach `%dopar%`- |
-
| 509 | -- |
- #' @importFrom foreach foreach- |
-
| 510 | -- |
- #' @importFrom ggplot2 aes_string- |
-
| 511 | -- |
- #' @importFrom ggplot2 geom_line- |
-
| 512 | -- |
- #' @importFrom ggplot2 geom_ribbon- |
-
| 513 | -- |
- #' @importFrom ggplot2 ggplot- |
-
| 514 | -- |
- #' @importFrom ggplot2 theme_bw- |
-
| 515 | -- |
- #' @importFrom ggplot2 scale_fill_manual- |
-
| 516 | -- |
- #' @importFrom ggplot2 xlab- |
-
| 517 | -- |
- #' @importFrom ggplot2 ylab- |
-
| 518 | -- |
- #' @importFrom magrittr "%>%"- |
-
| 519 | -- |
- #' @importFrom network network.size- |
-
| 520 | -- |
- #' @importFrom parallel makeCluster- |
-
| 521 | -- |
- #' @importFrom parallel stopCluster- |
-
| 522 | -- |
- #' @importFrom parallel clusterExport- |
-
| 523 | -- |
- #' @importFrom scales muted- |
-
| 524 | -- |
- #' @importFrom stats sd- |
-
| 525 | -- |
- #' @importFrom stats na.omit- |
-
| 526 | -- |
- #' @importFrom utils setTxtProgressBar- |
-
| 527 | -- |
- #' @importFrom utils txtProgressBar- |
-
| 528 | -- |
- #' @author Derek Corcoran <derek.corcoran.barrios@gmail.com>- |
-
| 529 | -- |
- #' @author M. Isidora Ávila-Thieme <msavila@uc.cl>- |
-
| 530 | -- |
- #' @author Erik Kusch <erik.kusch@bio.au.dk>- |
-
| 531 | -- |
- #' @export- |
-
| 532 | -- |
- RandomExtinctions <- function(Network, nsim = 10,- |
-
| 533 | -- |
- Record = FALSE, plot = FALSE,- |
-
| 534 | -- |
- SimNum = NULL,- |
-
| 535 | -- |
- NetworkType = "Trophic", clust.method = "cluster_infomap",- |
-
| 536 | -- |
- parallel = FALSE, ncores,- |
-
| 537 | -- |
- IS = 0,- |
-
| 538 | -- |
- Rewiring = FALSE, RewiringDist = NULL, RewiringProb = 0.5,- |
-
| 539 | -- |
- verbose = TRUE){
- |
-
| 540 | -! | -
- if(!NetworkType %in% c("Trophic", "Mutualistic")){stop("Please specify NetworkType as either 'Trophic' or 'Mutualistic'")}
- |
-
| 541 | -- |
- ## setting up objects- |
-
| 542 | -2x | -
- NumExt <- sd <- AccSecExt <- AccSecExt_95CI <- AccSecExt_mean <- Lower <- NULL- |
-
| 543 | -2x | -
- network <- .DataInit(x = Network)- |
-
| 544 | -2x | -
- if(is.null(SimNum)){
- |
-
| 545 | -2x | -
- SimNum <- network::network.size(network)- |
-
| 546 | -- |
- }- |
-
| 547 | -- | - - | -
| 548 | -- |
- ## simulations- |
-
| 549 | -2x | -
- if(verbose){ProgBar <- txtProgressBar(max = nsim, style = 3)}
- |
-
| 550 | -2x | -
- if(parallel){
- |
-
| 551 | -! | -
- cl <- makeCluster(ncores)- |
-
| 552 | -! | -
- registerDoParallel(cl)- |
-
| 553 | -! | -
- parallel::clusterExport(cl,- |
-
| 554 | -! | -
- varlist = c("network", "SimNum", "IS", "Rewiring", "RewiringDist", "RewiringProb"),
- |
-
| 555 | -! | -
- envir = environment()- |
-
| 556 | -- |
- )- |
-
| 557 | -! | -
- sims <- foreach(i=1:nsim, .packages = "NetworkExtinction")%dopar%{
- |
-
| 558 | -! | -
- sims <- try(ExtinctionOrder(Network = network, Order = sample(1:network::network.size(network), size = SimNum),- |
-
| 559 | -! | -
- IS = IS, NetworkType = NetworkType,- |
-
| 560 | -! | -
- Rewiring = Rewiring, RewiringDist = RewiringDist,- |
-
| 561 | -! | -
- verbose = FALSE, RewiringProb = RewiringProb), silent = TRUE)- |
-
| 562 | -! | -
- try({sims$simulation <- i}, silent = TRUE)
- |
-
| 563 | -! | -
- sims- |
-
| 564 | -- |
- }- |
-
| 565 | -! | -
- stopCluster(cl)- |
-
| 566 | -- |
- }else{
- |
-
| 567 | -2x | -
- sims <- list()- |
-
| 568 | -2x | -
- for(i in 1:nsim){
- |
-
| 569 | -120x | -
- sims[[i]] <- try(ExtinctionOrder(Network = network, Order = sample(1:network::network.size(network), size = SimNum),- |
-
| 570 | -120x | -
- IS = IS, NetworkType = NetworkType,- |
-
| 571 | -120x | -
- Rewiring = Rewiring, RewiringDist = RewiringDist,- |
-
| 572 | -120x | -
- verbose = FALSE, RewiringProb = RewiringProb), silent = TRUE)- |
-
| 573 | -120x | -
- try({sims[[i]]$simulation <- i}, silent = TRUE)
- |
-
| 574 | -120x | -
- if(verbose){setTxtProgressBar(ProgBar, i)}
- |
-
| 575 | -- |
- }- |
-
| 576 | -- |
- }- |
-
| 577 | -- | - - | -
| 578 | -- |
- ## extract objects- |
-
| 579 | -2x | -
- temps <- lapply(sims, "[[", 2)- |
-
| 580 | -2x | -
- sims <- lapply(sims, "[[", 1)- |
-
| 581 | -2x | -
- cond <- sapply(sims, function(x) "data.frame" %in% class(x))- |
-
| 582 | -2x | -
- cond <- c(1:length(cond))[cond]- |
-
| 583 | -2x | -
- sims <- sims[cond]- |
-
| 584 | -2x | -
- sims <- do.call(rbind, sims)- |
-
| 585 | -2x | -
- if(Record == TRUE){
- |
-
| 586 | -! | -
- FullSims <- sims- |
-
| 587 | -- |
- }- |
-
| 588 | -- | - - | -
| 589 | -2x | -
- sims <- sims[!is.na(sims$SecExt), ] %>% dplyr::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))- |
-
| 590 | -- | - - | -
| 591 | -- |
- ## plot output- |
-
| 592 | -2x | -
- if(plot == TRUE){
- |
-
| 593 | -! | -
- 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()
- |
-
| 594 | -! | -
- print(g)- |
-
| 595 | -- |
- }- |
-
| 596 | -- | - - | -
| 597 | -- |
- ## object output- |
-
| 598 | -2x | -
- if(Record == T & plot == T){
- |
-
| 599 | -! | -
- return(list(sims = sims, graph = g, FullSims = FullSims, nets = temps))- |
-
| 600 | -2x | -
- }else if(Record == F & plot == T){
- |
-
| 601 | -! | -
- return(list(sims = sims, graph = g, nets = temps))- |
-
| 602 | -2x | -
- }else if(Record == F & plot == F){
- |
-
| 603 | -2x | -
- return(list(sims = sims, nets = temps))- |
-
| 604 | -! | -
- }else if(Record == T & plot == F){
- |
-
| 605 | -! | -
- return(list(sims = sims, FullSims = FullSims, nets= temps))- |
-
| 606 | -- |
- }- |
-
| 607 | -- |
- }- |
-
| 608 | -- | - - | -
| 609 | -- |
- #' Comparison of Null hypothesis with other extinction histories- |
-
| 610 | -- |
- #'- |
-
| 611 | -- |
- #' It compares an object generated either by the Mostconnected or ExtinctionOrder functions- |
-
| 612 | -- |
- #' with a null hypothesis generated by the RandomExtinctions function it is important that- |
-
| 613 | -- |
- #' RandomExtinctions is in plot = T.- |
-
| 614 | -- |
- #'- |
-
| 615 | -- |
- #' @param Nullmodel an object generated by the RandomExtinctions- |
-
| 616 | -- |
- #' @param Hypothesis Extinction history generated by the Mostconnected or ExtinctionOrder- |
-
| 617 | -- |
- #' fuction- |
-
| 618 | -- |
- #' @return a plot comparing the expected value of secondary extinctions originated at random- |
-
| 619 | -- |
- #' with the observed extinction history.- |
-
| 620 | -- |
- #'- |
-
| 621 | -- |
- #' @examples- |
-
| 622 | -- |
- #' \dontrun{
- |
-
| 623 | -- |
- #' data("Less_Connected")
- |
-
| 624 | -- |
- #' History <- SimulateExtinctions(Network = Less_Connected, Method = "Mostconnected")- |
-
| 625 | -- |
- #' NullHyp <- RandomExtinctions(Network = Less_Connected, nsim = 100)- |
-
| 626 | -- |
- #' CompareExtinctions(Nullmodel = NullHyp, Hypothesis = History)- |
-
| 627 | -- |
- #' }- |
-
| 628 | -- |
- #' @importFrom broom tidy- |
-
| 629 | -- |
- #' @importFrom ggplot2 aes- |
-
| 630 | -- |
- #' @importFrom ggplot2 geom_line- |
-
| 631 | -- |
- #' @importFrom ggplot2 geom_point- |
-
| 632 | -- |
- #' @importFrom ggplot2 geom_ribbon- |
-
| 633 | -- |
- #' @importFrom ggplot2 scale_color_manual- |
-
| 634 | -- |
- #' @importFrom ggplot2 theme_bw- |
-
| 635 | -- |
- #' @importFrom scales muted- |
-
| 636 | -- |
- #' @author Derek Corcoran <derek.corcoran.barrios@gmail.com>- |
-
| 637 | -- |
- #' @author M. Isidora Ávila-Thieme <msavila@uc.cl>- |
-
| 638 | -- |
- #' @export- |
-
| 639 | -- | - - | -
| 640 | -- |
- CompareExtinctions <- function(Nullmodel, Hypothesis){
- |
-
| 641 | -1x | -
- if(class(Hypothesis$sims)[2] == "SimulateExt"){
- |
-
| 642 | -1x | -
- NumExt <- sd <- AccSecExt <- AccSecExt_mean <-NULL- |
-
| 643 | -1x | -
- if(class(Nullmodel$sims)[1] == "list"){
- |
-
| 644 | -! | -
- g <- Nullmodel$graph + geom_line(aes(color = "blue"))- |
-
| 645 | -! | -
- 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"))
- |
-
| 646 | -- |
- } else {
- |
-
| 647 | -1x | -
- g <- ggplot(Nullmodel$sims, aes(x = NumExt, y = AccSecExt_mean)) + geom_ribbon(aes_string(ymin = "Lower", ymax = "Upper"), fill = scales::muted("red")) + geom_line(color = "blue") + ylab("Acc. Secondary extinctions") + xlab("Primary extinctions") + theme_bw()
- |
-
| 648 | -1x | -
- 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"))
- |
-
| 649 | -1x | -
- g- |
-
| 650 | -- |
- }- |
-
| 651 | -- | - - | -
| 652 | -1x | -
- g- |
-
| 653 | -- | - - | -
| 654 | -1x | -
- return(g)- |
-
| 655 | -- |
- }- |
-
| 656 | -! | -
- if(class(Hypothesis$sims)[2] %in% c("Mostconnected", "ExtinctionOrder")){
- |
-
| 657 | -! | -
- NumExt <- sd <- AccSecExt <- AccSecExt_mean <-NULL- |
-
| 658 | -! | -
- g <- Nullmodel$graph + geom_line(aes(color = "blue"))- |
-
| 659 | -! | -
- 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"))
- |
-
| 660 | -! | -
- g- |
-
| 661 | -! | -
- return(g)- |
-
| 662 | -- |
- }- |
-
| 663 | -- |
- else{
- |
-
| 664 | -! | -
- message("Hypothesis not of class Mostconnected or ExtinctionOrder")
- |
-
| 665 | -- |
- }- |
-
| 666 | -- |
- }- |
-
| 1 | -- |
- #' Degree distribution of the network- |
-
| 2 | -- |
- #'- |
-
| 3 | -- |
- #' This function calculates the degree distribution of the network. First it- |
-
| 4 | -- |
- #' fits exponential, power law and truncated power law distribution models,- |
-
| 5 | -- |
- #' and calculates the AIC values to select the best fit, and finally it plots- |
-
| 6 | -- |
- #' the degree distribution in a log log scale showing the three fitted models- |
-
| 7 | -- |
- #' mentioned above against the observed distribution.- |
-
| 8 | -- |
- #'- |
-
| 9 | -- |
- #'- |
-
| 10 | -- |
- #' @param Network a trophic network of class network- |
-
| 11 | -- |
- #' @param scale a character stating if the graph is on a log-log scale- |
-
| 12 | -- |
- #' ("LogLog") or arithmetic scale ("arithmetic"), defaults to arithmetic
- |
-
| 13 | -- |
- #' @return exports three principal results:- |
-
| 14 | -- |
- #' 1. A list with network degree distribution values and with the value of each fit model- |
-
| 15 | -- |
- #' 2. A list with each model results and AIC of the distribution models- |
-
| 16 | -- |
- #' 3. A Ghraph of the degree distribution with the models adjust- |
-
| 17 | -- |
- #' In DDvalues, k represent the degree of the network and cumulative- |
-
| 18 | -- |
- #' the probability that each specie could be have this degree (pk).- |
-
| 19 | -- |
- #' Observation: In the graph, the zero values are not represented but this result are incorporate in the DF result- |
-
| 20 | -- |
- #'- |
-
| 21 | -- |
- #'@examples- |
-
| 22 | -- |
- #'library(NetworkExtinction)- |
-
| 23 | -- |
- #'data("chilean_intertidal")
- |
-
| 24 | -- |
- #'DegreeDistribution(chilean_intertidal)- |
-
| 25 | -- |
- #'- |
-
| 26 | -- |
- #'@importFrom sna degree- |
-
| 27 | -- |
- #'@importFrom stats nls- |
-
| 28 | -- |
- #'@importFrom broom augment- |
-
| 29 | -- |
- #'@importFrom broom glance- |
-
| 30 | -- |
- #'@importFrom dplyr arrange- |
-
| 31 | -- |
- #'@importFrom dplyr bind_rows- |
-
| 32 | -- |
- #'@importFrom dplyr case_when- |
-
| 33 | -- |
- #'@importFrom dplyr filter- |
-
| 34 | -- |
- #'@importFrom dplyr full_join- |
-
| 35 | -- |
- #'@importFrom dplyr group_split- |
-
| 36 | -- |
- #'@importFrom dplyr mutate- |
-
| 37 | -- |
- #'@importFrom dplyr select- |
-
| 38 | -- |
- #'@importFrom tidyr gather- |
-
| 39 | -- |
- #'@importFrom magrittr "%>%"- |
-
| 40 | -- |
- #'@importFrom ggplot2 aes_string- |
-
| 41 | -- |
- #'@importFrom ggplot2 geom_line- |
-
| 42 | -- |
- #'@importFrom ggplot2 geom_point- |
-
| 43 | -- |
- #'@importFrom ggplot2 scale_x_log10- |
-
| 44 | -- |
- #'@importFrom ggplot2 scale_y_log10- |
-
| 45 | -- |
- #'@importFrom ggplot2 theme_bw- |
-
| 46 | -- |
- #'@importFrom ggplot2 ylim- |
-
| 47 | -- |
- #'@importFrom MASS fitdistr- |
-
| 48 | -- |
- #'@importFrom purrr map- |
-
| 49 | -- |
- #'@importFrom purrr reduce- |
-
| 50 | -- |
- #'@importFrom stats ks.test- |
-
| 51 | -- |
- #'@importFrom stats glm- |
-
| 52 | -- |
- #'@importFrom stats logLik- |
-
| 53 | -- |
- #'@importFrom stats predict- |
-
| 54 | -- |
- #' @author Derek Corcoran <derek.corcoran.barrios@gmail.com>- |
-
| 55 | -- |
- #' @author M.Isidora Avila Thieme <msavila@uc.cl>- |
-
| 56 | -- |
- #' @export- |
-
| 57 | -- | - - | -
| 58 | -- | - - | -
| 59 | -- |
- DegreeDistribution <- function(Network, scale = "arithmetic"){
- |
-
| 60 | -1x | -
- AIC <- Cumulative <- Exp <- fit <- model <- LogPower <- logLik <- BIC <- Power <- Normal.Resid <- LogExp <- family <- AICcNorm <- NULL- |
-
| 61 | -1x | -
- Network <- .DataInit(Network)- |
-
| 62 | -1x | -
- totaldegree<- degree(Network)- |
-
| 63 | -1x | -
- K <- 0:max(totaldegree)- |
-
| 64 | -1x | -
- For.Graph<- data.frame(K = K, Cumulative = NA)- |
-
| 65 | -1x | -
- for(i in 1:length(K)){
- |
-
| 66 | -68x | -
- For.Graph$Cumulative[i] <- sum(totaldegree>K[i])/length(totaldegree)- |
-
| 67 | -- |
- }- |
-
| 68 | -1x | -
- For.Graph <- For.Graph %>% mutate(LogK = log(K), LogCum = log(Cumulative))- |
-
| 69 | -- | - - | -
| 70 | -- | - - | -
| 71 | -- |
- #exponential model nls- |
-
| 72 | -1x | -
- exp.model <- nls(Cumulative~exp(K*lambda+ c),start= list(lambda=0.1, c = 0), data = For.Graph)- |
-
| 73 | -1x | -
- For.Graph$Exp <- predict(exp.model)- |
-
| 74 | -1x | -
- Summs.exp <- glance(exp.model)- |
-
| 75 | -1x | -
- Summs.exp$model <- "Exp"- |
-
| 76 | -1x | -
- Summs.exp$Normal.Resid <- ifelse(tidy(ks.test(augment(exp.model)$.resid,y='pnorm',alternative='two.sided'))$p.value < 0.05, "No", "Yes")- |
-
| 77 | -1x | -
- Summs.exp$family <- "Exponential"- |
-
| 78 | -1x | -
- Summs.exp$AICcNorm <- logLik(MASS::fitdistr(augment(exp.model)$.resid, "normal"))[1]- |
-
| 79 | -1x | -
- Summs.exp$AICcNorm <- (4 - 2*Summs.exp$AICcNorm) + (12/(nrow(augment(exp.model)) - 1))- |
-
| 80 | -1x | -
- Params.exp <- tidy(exp.model)- |
-
| 81 | -1x | -
- Params.exp$model <- "Exp"- |
-
| 82 | -- |
- #exponential model Log- |
-
| 83 | -- | - - | -
| 84 | -1x | -
- power <- filter(For.Graph, K != 0 & Cumulative != 0)- |
-
| 85 | -1x | -
- logexp.model <- glm(LogCum ~ K, data = power)- |
-
| 86 | -1x | -
- power$LogExp <- exp(predict(logexp.model))- |
-
| 87 | -1x | -
- Summs.logexp <- glance(logexp.model)- |
-
| 88 | -1x | -
- Summs.logexp$model <- "LogExp"- |
-
| 89 | -1x | -
- Summs.logexp$Normal.Resid <- ifelse(tidy(ks.test(augment(logexp.model)$.resid,y='pnorm',alternative='two.sided'))$p.value < 0.05, "No", "Yes")- |
-
| 90 | -1x | -
- Summs.logexp$family <- "Exponential"- |
-
| 91 | -1x | -
- Summs.logexp$AICcNorm <- logLik(MASS::fitdistr(augment(logexp.model)$.resid, "normal"))[1]- |
-
| 92 | -1x | -
- Summs.logexp$AICcNorm <- (4 - 2*Summs.logexp$AICcNorm) + (12/(nrow(augment(logexp.model)) - 1))- |
-
| 93 | -1x | -
- Params.logexp <- tidy(logexp.model)- |
-
| 94 | -1x | -
- Params.logexp$model <- "LogExp"- |
-
| 95 | -- | - - | -
| 96 | -- |
- #logpowerlaw- |
-
| 97 | -- | - - | -
| 98 | -1x | -
- logpower.model <- glm(LogCum ~ I(log(K)), data = power)- |
-
| 99 | -1x | -
- power$LogPower <- exp(predict(logpower.model))- |
-
| 100 | -1x | -
- For.Graph <- full_join(For.Graph, power)- |
-
| 101 | -1x | -
- Summs.logpower <- glance(logpower.model)- |
-
| 102 | -1x | -
- Summs.logpower$model <- "LogPower"- |
-
| 103 | -1x | -
- Summs.logpower$Normal.Resid <- ifelse(tidy(ks.test(augment(logpower.model)$.resid,y='pnorm',alternative='two.sided'))$p.value < 0.05, "No", "Yes")- |
-
| 104 | -1x | -
- Summs.logpower$family <- "PowerLaw"- |
-
| 105 | -1x | -
- Summs.logpower$AICcNorm <- logLik(MASS::fitdistr(augment(logpower.model)$.resid, "normal"))[1]- |
-
| 106 | -1x | -
- Summs.logpower$AICcNorm <- (4 - 2*Summs.logpower$AICcNorm) + (12/(nrow(augment(logpower.model)) - 1))- |
-
| 107 | -1x | -
- Params.logpower <- tidy(logpower.model)- |
-
| 108 | -1x | -
- Params.logpower$model <- "LogPower"- |
-
| 109 | -- | - - | -
| 110 | -- |
- #powerlaw- |
-
| 111 | -- | - - | -
| 112 | -1x | -
- powerlaw.model <- nls(Cumulative~a*K^y, start= list(y=0, a = 1), data = power)- |
-
| 113 | -1x | -
- power$Power <- predict(powerlaw.model)- |
-
| 114 | -1x | -
- For.Graph <- full_join(For.Graph, power)- |
-
| 115 | -1x | -
- Summs.power <- glance(powerlaw.model)- |
-
| 116 | -1x | -
- Summs.power$model <- "Power"- |
-
| 117 | -1x | -
- Summs.power$Normal.Resid <- ifelse(tidy(ks.test(augment(powerlaw.model)$.resid,y='pnorm',alternative='two.sided'))$p.value < 0.05, "No", "Yes")- |
-
| 118 | -1x | -
- Summs.power$family <- "PowerLaw"- |
-
| 119 | -1x | -
- Summs.power$AICcNorm <- logLik(MASS::fitdistr(augment(powerlaw.model)$.resid, "normal"))[1]- |
-
| 120 | -1x | -
- Summs.power$AICcNorm <- (4 - 2*Summs.power$AICcNorm) + (12/(nrow(augment(powerlaw.model)) - 1))- |
-
| 121 | -1x | -
- Params.power <- tidy(powerlaw.model)- |
-
| 122 | -1x | -
- Params.power$model <- "Power"- |
-
| 123 | -- | - - | -
| 124 | -- |
- #all together- |
-
| 125 | -1x | -
- Summs <- full_join(Summs.exp, Summs.power)- |
-
| 126 | -1x | -
- Summs <- full_join(Summs, Summs.logexp)- |
-
| 127 | -1x | -
- Summs <- full_join(Summs, Summs.logpower) %>% select(logLik, AIC, BIC, model, Normal.Resid, family, AICcNorm)- |
-
| 128 | -1x | -
- Summs <- arrange(Summs, Normal.Resid, AIC) %>% dplyr::select(logLik, AIC, BIC, model, Normal.Resid, family)- |
-
| 129 | -1x | -
- params <- bind_rows(Params.logpower, Params.power, Params.logexp, Params.exp) %>%- |
-
| 130 | -1x | -
- dplyr::filter(model %in% Summs$model) %>%- |
-
| 131 | -1x | -
- mutate(term = case_when(term == "y" ~ "Beta",- |
-
| 132 | -1x | -
- term == "a" ~ "c",- |
-
| 133 | -1x | -
- term == "(Intercept)" ~ "c",- |
-
| 134 | -1x | -
- term == "I(log(K))" ~ "Beta",- |
-
| 135 | -1x | -
- term == "lambda" ~ "Lambda",- |
-
| 136 | -1x | -
- term == "K" ~ "Lambda",- |
-
| 137 | -1x | -
- TRUE ~ term))- |
-
| 138 | -1x | -
- DF2 <- For.Graph %>% filter(K != 0 & Cumulative != 0) %>% gather(key = model, value = fit, Exp, Power, LogExp, LogPower) %>% dplyr::filter(model %in% Summs$model)- |
-
| 139 | -- | - - | -
| 140 | -1x | -
- g <- ggplot(DF2, aes_string(x = "K", y = "Cumulative")) + geom_line() + geom_point()+ theme_bw() + geom_line(aes_string(y ="fit", color = "model")) + ylim(c(0,1))- |
-
| 141 | -- | - - | -
| 142 | -1x | -
- if(scale == "LogLog"){
- |
-
| 143 | -! | -
- g <- g + scale_x_log10() + scale_y_log10(breaks=c(0, .001,.01,1))- |
-
| 144 | -- |
- }- |
-
| 145 | -- | - - | -
| 146 | -1x | -
- g- |
-
| 147 | -- | - - | -
| 148 | -1x | -
- return(list(DDvalues = For.Graph, models = Summs, graph = g, params = params))- |
-
| 149 | -- |
- }- |
-
| 150 | -- | - - | -
| 151 | -- | - - | -
| 1 | -- | - - | -
| 2 | -- |
- #' @param x a network representation as a an adyacency matrix, edgelist,- |
-
| 3 | -- |
- #' or a network object- |
-
| 4 | -- |
- #' @importFrom network as.network- |
-
| 5 | -- |
- #' @importFrom methods is- |
-
| 6 | -- |
- #' @noRd- |
-
| 7 | -- | - - | -
| 8 | -- |
- .DataInit <- function(x){
- |
-
| 9 | -140x | -
- if(is(x) == "network"){
- |
-
| 10 | -140x | -
- x- |
-
| 11 | -- |
- }- |
-
| 12 | -140x | -
- if(is(x) == "matrix"){
- |
-
| 13 | -! | -
- x <- as.network(x, loops = TRUE)- |
-
| 14 | -- |
- }- |
-
| 15 | -140x | -
- return(x)- |
-
| 16 | -- |
- }- |
-
| 1 | -- |
- #' Plots the extinctions history of a network- |
-
| 2 | -- |
- #'- |
-
| 3 | -- |
- #' It takes a NetworkTopology class object and plots the network index- |
-
| 4 | -- |
- #' after every extinction- |
-
| 5 | -- |
- #'- |
-
| 6 | -- |
- #' @param History a NetworkTopology object obtained from the Mostconnected function- |
-
| 7 | -- |
- #' or the ExtinctionOrder function- |
-
| 8 | -- |
- #' @param Variable the variable of the NetworkTopology object that you want as a y variable- |
-
| 9 | -- |
- #' @return A plot of number of extinctions in the x axis vs the choosen variable in the Y axis- |
-
| 10 | -- |
- #' @examples- |
-
| 11 | -- |
- #' # If you don't specify the y variable it will plot the secondary extinctions- |
-
| 12 | -- |
- #' # by default- |
-
| 13 | -- |
- #' data("net")
- |
-
| 14 | -- |
- #' history <- SimulateExtinctions(Network = net, Method = "Mostconnected")- |
-
| 15 | -- |
- #' ExtinctionPlot(History = history$sims)- |
-
| 16 | -- |
- #' # You can also specify the variable to be ploted in the y axis- |
-
| 17 | -- |
- #' ExtinctionPlot(History = history$sims, Variable = "Link_density")- |
-
| 18 | -- |
- #' @importFrom ggplot2 aes_string- |
-
| 19 | -- |
- #' @importFrom ggplot2 geom_line- |
-
| 20 | -- |
- #' @importFrom ggplot2 ggplot- |
-
| 21 | -- |
- #' @importFrom ggplot2 theme_bw- |
-
| 22 | -- |
- #' @author Derek Corcoran <derek.corcoran.barrios@gmail.com>- |
-
| 23 | -- |
- #' @author M. Isidora Ávila-Thieme <msavila@uc.cl>- |
-
| 24 | -- |
- #' @seealso [NetworkExtintion::ExtinctionOrder()]- |
-
| 25 | -- |
- #' @export- |
-
| 26 | -- | - - | -
| 27 | -- |
- ExtinctionPlot <- function(History, Variable = "AccSecExt"){
- |
-
| 28 | -1x | -
- History$X <- 1:nrow(History)- |
-
| 29 | -1x | -
- ggplot(History, aes_string(x = "X", y = Variable)) + geom_line() + theme_bw() + ylab(Variable) + xlab("Primary extinctions")
- |
-
| 30 | -- |
- }- |
-
-
+
The goal of NetworkExtinction is to Simulate the extinction of species in the food web and to analyze its cascading effects, as described in Dunne et al. (2002) \doi{10.1073/pnas.192407699}
diff --git a/README.md b/README.md
index a75b609..faccf9c 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,8 @@
+[](https://zenodo.org/badge/latestdoi/116978043)
+
# NetworkExtinction
# pkgdown
@@ -63,7 +65,16 @@ history <- SimulateExtinctions(Network = net, Method = "Mostconnected")
ExtinctionPlot(History = history[[1]], Variable = "AccSecExt")
```
+
++Figure 3. The graph shows the number of accumulated secondary +extinctions that occur when removing species from the most to the least +connected species +
+ +
++Figure 4. The graph shows the number of accumulated secondary +extinctions that occur when removing species in a custom order. +
+ +
++Figure 5. The resulting graph of the RandomExtinctions function +
+ +
++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 +
+ +
++Figure 7. Example of the use of the ExtinctionPlot function showing the +accumulated secondary extinctions against number of extinctions +
+ +
++Figure 8. Another example of the use of the ExtinctionPlot function +showing the number of links per species against number of extinctions +
+ +
++Figure 9: Fitted vs observed values of the degree distribution. The +black line and points show the observed values, the red, green and blue +lines show the fitted values for the Exponential, power law and trucated +distribution, respectively +
+ +