In this script I build random network with varying link distribution. Two things are being varied: the type of distribution (ie normal, poisson, uniform) and the parameters of these distribution (ie mean, upper range, lambda …). I built two-level network like plant-herbivores with 60 prey species and 200 predator species. I varied the threshold from which a secondary extinction occure, ie a predator species is removed from the network if 50% of its resources has gone extinct. I simulate random prey species extinction.
Below is the script:
# libraries needed
library(plyr)
library(dplyr)
##
## Attaching package: 'dplyr'
##
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(truncnorm) #for the truncated normal distribution
# these are the function doing the job a function to generate a random
# network in this version only two levels
create_net <- function(n1 = 20, n2 = 50, distr, params) {
if (length(distr) != n2)
stop
if (distr == "rnorm")
link_d <- rnorm(n2, params[1], params[2])
if (distr == "runif")
link_d <- runif(n2, params[1], params[2])
if (distr == "rpois")
link_d <- rpois(n2, params[1])
if (distr == "rtruncnorm")
link_d <- rtruncnorm(n2, a = params[1], b = params[2], mean = params[3],
sd = params[4])
# can implement other distribution
link_d <- round(link_d, 0)
link_d[link_d <= 0] <- 1
link_d[link_d > n1] <- n1
# put ones in the network
ones <- lapply(link_d, function(x) sample(x = 1:n1, size = x, replace = FALSE))
zeros <- rep(0, n1)
mat <- sapply(ones, function(x) {
zeros[x] <- 1
return(zeros)
}, simplify = "array")
return(mat)
}
# next step the function that remove randomly rows and check for secondary
# extinction based on certain criteria
extinction <- function(network, thresh = 0.5) {
link_d <- colSums(network)
step <- data.frame(N1 = dim(network)[1]:0, N2 = 0, Thresh = thresh)
step[1, 2] <- dim(network)[2]
# as long as the prey species are more than 2 do the loop
while (dim(network)[1] > 2) {
# randomly remove one prey species
network <- network[-sample(x = 1:dim(network)[1], size = 1), ]
new_link <- colSums(network)
link_ratio <- new_link/link_d
# if the remaining links per predator species are below the threshold remove
# the predator species
step[step$N1 == dim(network)[1], 2] <- length(which(link_ratio > thresh))
network <- network[, -link_ratio <= thresh]
if (dim(network)[2] == 0)
break
}
return(step)
}
# wrap the extinction function over nsim simulation, also re-creating the
# networks
ext_sim <- function(network, thresh = 0.5, nsim = 10) {
li <- lapply(1:nsim, function(x) extinction(network, thresh = thresh))
li <- ldply(li, function(x) x)
li$Simu <- rep(1:nsim, each = dim(network)[1] + 1)
# get the mean and standard deviation across the simulation
li_d <- li %>% group_by(N1) %>% summarise(N2_m = mean(N2), N2_sd = sd(N2),
Thresh = unique(Thresh))
return(li_d)
}
# now another wrap to test a range of thresholds
ext_range <- function(network, thresh = seq(0.1, 1, 0.3), nsim = 10) {
li <- lapply(thresh, function(x) ext_sim(network, x, nsim))
li <- ldply(li, function(x) x)
return(li)
}
#compare some random normal network
net1<-create_net(n1=60,n2=200,distr="rnorm",params=c(30,5))
net2<-create_net(n1=60,n2=200,distr="rnorm",params=c(30,30))
net3<-create_net(n1=60,n2=200,distr="rnorm",params=c(2,5))
net4<-create_net(n1=60,n2=200,distr="rnorm",params=c(50,5))
net_li<-list(net1,net2,net3,net4)
net_res<-lapply(net_li,function(x) ext_range(x))
dim_res<-lapply(net_res,function(x) dim(x))
net_res<-ldply(net_res,function(x) x)
net_res$Type<-rep(c("Large_M_Low_SD","Large_M_Large_SD","Low_M_Low_SD","VLarge_M_Low_SD"),each=244)
ggplot(net_res,aes(x=N1,y=N2_m,color=factor(Thresh)))+geom_point()+geom_errorbar(aes(ymin=N2_m-N2_sd,ymax=N2_m+N2_sd))+facet_wrap(~Type)+labs(x="Plant Species Richness",y="Consumer Species Richness",title="Secondary extinction with varying threshold\nand different network creation from a normal distribution")
#same stuff for uniform distribution
net1<-create_net(n1=60,n2=200,distr="runif",params=c(1,60))
net2<-create_net(n1=60,n2=200,distr="runif",params=c(35,45))
net3<-create_net(n1=60,n2=200,distr="runif",params=c(1,15))
net4<-create_net(n1=60,n2=200,distr="runif",params=c(45,60))
net_li<-list(net1,net2,net3,net4)
net_res<-lapply(net_li,function(x) ext_range(x))
net_res<-ldply(net_res,function(x) x)
net_res$Type<-rep(c("Wide","Restricted_center","Restricted_left","Restricted_right"),each=244)
ggplot(net_res,aes(x=N1,y=N2_m,color=factor(Thresh)))+geom_point()+geom_errorbar(aes(ymin=N2_m-N2_sd,ymax=N2_m+N2_sd))+facet_wrap(~Type)+labs(x="Plant Species Richness",y="Consumer Species Richness",title="Secondary extinction with varying threshold\nand different network creation from a uniform distribution")
#same stuff for poisson network
net1<-create_net(n1=60,n2=200,distr="rpois",params=1)
net2<-create_net(n1=60,n2=200,distr="rpois",params=10)
net3<-create_net(n1=60,n2=200,distr="rpois",params=30)
net4<-create_net(n1=60,n2=200,distr="rpois",params=50)
net_li<-list(net1,net2,net3,net4)
net_res<-lapply(net_li,function(x) ext_range(x))
net_res<-ldply(net_res,function(x) x)
net_res$Type<-rep(c("Restricted","Low","Medium","High"),each=244)
ggplot(net_res,aes(x=N1,y=N2_m,color=factor(Thresh)))+geom_point()+geom_errorbar(aes(ymin=N2_m-N2_sd,ymax=N2_m+N2_sd))+facet_wrap(~Type)+labs(x="Plant Species Richness",y="Consumer Species Richness",title="Secondary extinction with varying threshold\nand different network creation from a poisson distribution")
#same stuff for truncated network
net1<-create_net(n1=60,n2=200,distr="rtruncnorm",params=c(1,60,30,10))
net2<-create_net(n1=60,n2=200,distr="rtruncnorm",params=c(1,10,3,4))
net3<-create_net(n1=60,n2=200,distr="rtruncnorm",params=c(40,60,55,10))
net4<-cbind(create_net(n1=60,n2=100,distr="rtruncnorm",params=c(1,10,3,4)),create_net(n1=60,n2=100,distr="rtruncnorm",params=c(40,60,55,10)))
net_li<-list(net1,net2,net3,net4)
net_res<-lapply(net_li,function(x) ext_range(x))
net_res<-ldply(net_res,function(x) x)
net_res$Type<-rep(c("Wide","Left","Right","Bimodal"),each=244)
ggplot(net_res,aes(x=N1,y=N2_m,color=factor(Thresh)))+geom_point()+geom_errorbar(aes(ymin=N2_m-N2_sd,ymax=N2_m+N2_sd))+facet_wrap(~Type)+labs(x="Plant Species Richness",y="Consumer Species Richness",title="Secondary extinction with varying threshold\nand different network creation from a truncated normal distribution")
Several points here: 1. Contrary to my expectations a uniform distribution of links to not lead to different patterns as the one observed with the normal distribution 2. Also surprising is the fact that specialist are preventing the collapse of the network. In networks with only generalist (high number of links, low variation around this high number) all predator species go extinct at the same level of primary extinction. As we are doing random primary extinction, all generalist will reach the threshold at similar point leading to the observed collapse. 3. On the opposite in networks made mostly of specialist we observe a linearization of the secondary extinction. This makes sense as every species is linked to different host under random host removal only a few of them will be affected at each host removal step.
So very exciting insight for me, I can see several next step from this: