Nothing
R/networklevel.R
In bipartite: Visualising Bipartite Networks and Calculating Some (Ecological) Indices
Defines functions
.networklevel
`networklevel`
`networklevel` <- function(web, index="ALLBUTDD", level="both", weighted=TRUE, ISAmethod="Bluethgen", SAmethod="Bluethgen", extinctmethod="r", nrep=100, CCfun=median, dist="horn", normalise=TRUE, empty.web=TRUE, logbase="e", intereven="prod", H2_integer=TRUE, fcweighted=TRUE, fcdist="euclidean", legacy=FALSE){
##
## web interaction matrix, with lower trophic level in rows, higher in columns
## legacy se to TRUE allows to run networklevel in its old form
#! JFedit: fddist and fdweighted replace throughout by fcdist and fcweighted
## !!!!! when adding a new index, make sure to also add it to the "exclude from grouplevel"-list !!!!!##
if(empty.web) {web <- empty(web)}
web.e <- empty(web) # emptied web for some indices
if (NROW(web) < 2 | NCOL(web) <2) warning("Web is really too small to calculate any reasonable index. You will get the values nonetheless, but I wouldn't put any faith in them!")
allindex <- c( #descriptive:
"number of species", "connectance", "web asymmetry",
#binary based and related:
"links per species", "number of compartments", "compartment diversity",
"cluster coefficient", "degree distribution", "mean number of shared partners",
"togetherness", "C score", "V ratio", "discrepancy", "nestedness", "NODF", "weighted nestedness",
#miscelleneous:
"ISA", "SA", "extinction slope", "robustness", "niche overlap",
#quantitative series:
"weighted cluster coefficient", "modularity","weighted NODF", "partner diversity", "generality", "vulnerability", "linkage density", "weighted connectance", "Fisher alpha", "interaction evenness", "Alatalo interaction evenness", "Shannon diversity", "functional complementarity", "H2" )
# GONE: "mean interaction diversity", "effective partners",
index <- unique(index) # enforces that each index name is used only once
wrong.name <- which(is.na(pmatch(index, c(allindex, "ALL", "ALLBUTDD", "info","quantitative", "binary", "topology", "networklevel"))))
if (length(wrong.name) > 0) stop("You selected an index that is not available: ", paste(index[wrong.name], collapse=", "))
# only if indices are not given explicitly:
if (length(index)==1 & !all(index %in% allindex)){
index <- switch(index,
"ALL" = allindex,
"ALLBUTDD" = allindex[-which(allindex=="degree distribution")],
"info" = c("number of species", "connectance", "web asymmetry", "links per species", "number of compartments"),
# logic: only rough information on the network's general structure
"quantitative" = c("weighted cluster coefficient", "modularity", "weighted nestedness", "weighted NODF", "functional complementarity", "partner diversity", "effective partners", "H2", "diversity","linkage density", "weighted connectance", "niche overlap"), #"mean interaction diversity",
# logic: the "quantitative series"
"binary" = c("connectance", "links per species", "nestedness", "mean number of partners","cluster coefficient", "C-score", "Fisher alpha"),
# logic: metrics for binary networks
"topology" = c("connectance", "cluster coefficient", "degree distribution", "togetherness", "nestedness", "NODF"),
# logic: more abstract, topological metrics for binary networks
"networklevel" = c("connectance", "web asymmetry", "links per species", "number of compartments", "compartment diversity", "cluster coefficient", "modularity", "nestedness", "NODF", "weighted NODF", "ISA", "SA", "linkage density", "Fisher alpha", "diversity", "interaction evenness", "Alatalo interaction evenness", "H2"),
# only the truly networky indices
stop("Your index is not recognised! Typo? Check help for options!", call.=FALSE) #default for all non-matches
)
}
if (legacy == FALSE){
out <- list()
# first: all the true networklevel indices
#--------------------
if ("connectance" %in% index){
# connectance: "the fraction of all possible links that are realized in a network", (p. 12917 in Dunne et al. 2002)
suppressWarnings(out$connectance <- sum(web>0)/prod(dim(web)))
}
#--------------------
if ("web asymmetry" %in% index) out$"web asymmetry" <- (NCOL(web)-NROW(web))/sum(dim(web)) # web asymmetry (Bluethgen et al. 2007, Fig. S2)
#--------------------
if ("links per species" %in% index){
L <- sum(web>0)/sum(dim(web))
out$"links per species"=L
}
#--------------------
if (any(c("number of compartments", "compartment diversity") %in% index)){
CD <- function(co){
if (co$n.compart>1){
no <- NA
for (i in 1:co$n.compart){
comp <- which(abs(co$cweb)==i, arr.ind=TRUE) # who is in a compartment?
no[i] <- length(unique(comp[,1])) + length(unique(comp[,2])) # how many species
}
no <- no/sum(dim(web)) # standardise for number of species in the web
CD <- exp(-sum(no*log(no)))
} else {CD <- NA} #; warning("only one compartment")}
CD
}
comps <- try(compart(web.e), silent=TRUE)
if (inherits(comps, "try-error")) {
ncompart <- compdiv <- NA
} else {
ncompart <- comps$n.compart
compdiv <- CD(comps)
}
if ("number of compartments" %in% index) out$"number of compartments" <- as.integer(ncompart)
if ("compartment diversity" %in% index) out$"compartment diversity" <- compdiv
}
#------------------------
if ("cluster coefficient" %in% index){
cluster.coef <- function(web, full=FALSE, FUN=mean){
# calculate cluster coefficient
# web a bipartite web
# full logical; return cluster coefficients for each species in the network?
# FUN give a function to summarise individual cluster coefficients, defaults to 'mean'.
# The concept was introduced by Watts & Strogatz (1998) and is described under in Wikipedia under http://en.wikipedia.org/w/index.php?title=Clustering_coefficient
# Its main idea was to help identifying small-world networks, which should
# have high cluster coefficients but low path lengths.
# Ci of species i is simply the number of realised links devided by the number of possible links. At the network level, the CC for a network is the average CC of its members. This is a little bit fishy, since the CCs are log-normal
# distributed (in pollinations networks at least). Therefore with FUN also
# other summary measures can be employed.
# Because within 'bipartite' we look at 2-mode networks, mean C is:
# C_lowertrophiclevel = C_highertrophiclevel = C_entirenetwork.
#
# Literature: Watts DJ, Strogatz S (1998) Collective dynmics of 'small-world' networks. Nature 393:440-442
# notice that this CC differs from the binary version of Tore Opsahl's. This here is the same as the one-mode cluster coefficient, Tore's is the "corrected" two-mode!
web <- as.matrix(web)
Ci.high <- colSums((web>0))/nrow(web)
Ci.low <- rowSums((web>0))/ncol(web)
CC <- FUN(Ci.high)
if (full) out <- list("cluster coefficient"=CC, "CC values higher"=Ci.high,
"CC values lower"=Ci.low) else out <- c("cluster coefficient"=CC)
out
}
out$"cluster coefficient"=as.numeric(cluster.coef(web, FUN=CCfun, full=FALSE))
}
#-------------------
if ("modularity" %in% index){
Q <- try(computeModules(web)@likelihood)
out$"modularity Q" <- ifelse(inherits(Q, "try-error"), NA, Q)
}
#-------------------
if ("nestedness" %in% index){
nest <- try(nestedtemp(web)$statistic, silent=TRUE)
out$nestedness <- ifelse(inherits(nest, "try-error"), NA, nest)
# a fast implementation of nestedness by Jari Oksanen
#old: nestedness(web, null.models=FALSE)$temperature
}
#-------------------
if ("NODF" %in% index){
NODF <- try(unname(nestednodf(web, order=TRUE, weighted=FALSE)$statistic[3]), silent=TRUE)
out$"NODF" <- if (inherits(NODF, "try-error")) NA else NODF
}
#-------------------
if ("weighted nestedness" %in% index){
wine.res <- try(wine(web.e, nreps=nrep)$wine, silent=TRUE)
out$"weighted nestedness" <- if (!inherits(wine.res, "try-error")) {wine.res} else {NA}
}
#-------------------
if ("weighted NODF" %in% index){
wNODF <- try(unname(nestednodf(web, order=TRUE, weighted=TRUE)$statistic[3]), silent=TRUE)
out$"weighted NODF" <- if (inherits(wNODF, "try-error")) NA else wNODF
}
#------------------
if (any(c("ISA", "interaction strength asymmetry", "dependence asymmetry") %in% index)){
# Dependence asymmetry (Bascompte et al. 2006; Bluethgen et al. 2007, Fig. S2)
depL <- web.e/matrix(rowSums(web.e), nrow=NROW(web.e), ncol=NCOL(web.e), byrow=FALSE)
depH <- web.e/matrix(colSums(web.e), nrow=NROW(web.e), ncol=NCOL(web.e), byrow=TRUE)
if (ISAmethod=="Bascompte" & "ISA" %in% index) {
#depMax <- depL
#greaterindex <- depL < depH
#depMax[greaterindex] <- depH[greaterindex]
out$"dependence asymmetry"=mean(abs(depL-depH)/pmax(depL, depH), na.rm=TRUE)
}
if (ISAmethod=="Bluethgen" & "ISA" %in% index) {
web2 <- web
# delete cells for species encountered only once (Bluethgen, pers. comm.):
web2[, which(colSums(web)==1)] <- 0
web2[which(rowSums(web)==1), ] <- 0
rowsummat <- matrix(rowSums(web2), nrow=NROW(web2), ncol=NCOL(web2), byrow=FALSE)
colsummat <- matrix(colSums(web2), nrow=NROW(web2), ncol=NCOL(web2), byrow=TRUE)
depL <- web2/rowsummat
depH <- web2/colsummat
depL[depL<=0] <- NA
depH[depH<=0] <- NA
# now we need a correction to account for the fact that links with few (e.g. 2) observations will have a minimum depL of 1/2: all on one species: depL=1, one on each of two: depL=0.5
depLprime <- (depL - 1/rowsummat)/(1 - 1/rowsummat)
# assumes depLmin = 1/web2 and depLmax=1
depHprime <- (depH - 1/colsummat)/(1 - 1/colsummat)
out$"interaction strength asymmetry"=mean(as.matrix(depHprime-depLprime), na.rm=TRUE) #ranges from -1 to 1 /sum(depLprime, depHprime, na.rm=TRUE)
}
}
#------------------
if ("SA" %in% index){
# Specialisation asymmetry (Bluethgen et al. 2007, Fig. S2)
# 2 options for calculating the "mean" SA:
# either as Bluethgen et al: average weighted by number of interactions in the
# cell or as mean of logarithms (since the dependencies follow a lognormal
# distribution)
di <- dfun(web)$dprime # plants
dj <- dfun(t(web))$dprime # pollinators
if (SAmethod=="log"){
lgmeani <- mean(log(di[di>0])); lgmeanj <- mean(log(dj[dj>0]))
SA <- (lgmeanj-lgmeani)/sum(lgmeani, lgmeanj)
# ij-sequence changed because log changes sequence, too
}
if (SAmethod=="Bluethgen"){
wmeani <- sum(di*rowSums(web.e))/sum(web.e)
wmeanj <- sum(dj*colSums(web.e))/sum(web.e)
SA <- (wmeanj-wmeani)/sum(wmeani, wmeanj)
# positive values indicate more specialisation in the higher trophic level
}
out$"specialisation asymmetry" <- SA
}
#------------------
if (any(c("linkage density","weighted connectance") %in% index)){
# for formula see Tylianakis et al. (2006), supplement.
# N refers to prey, P to predators
preytot.mat <- matrix(rep(colSums(web), NROW(web)), NROW(web), byrow=TRUE)
preyprop.mat <- web/preytot.mat # = b_ik/b_.k in the first formula
#H_Nk is the diversity index of inflow (diversity of flower visits for each pollinator)
predtot.mat <- matrix(rep(rowSums(web), NCOL(web)), NROW(web), byrow=FALSE)
predprop.mat <- web/predtot.mat # = b_kj/b_.k in the second formula
if (logbase==2 | logbase=="2"){
H_Nk <- apply(preyprop.mat, 2, function(x) -sum(x*log2(x), na.rm=TRUE))
#H_Pk is the diversity index of pollinators for each plant species
H_Pk <- apply(predprop.mat, 1, function(x) -sum(x*log2(x), na.rm=TRUE))
# next, we need the reciprocals of this
# note that the ifelse is only needed if the web contains prey that is
# not eaten or predators that don't eat ...
n_Nk <- ifelse(colSums(web)!=0, 2^H_Nk, 0)
n_Pk <- ifelse(rowSums(web)!=0, 2^H_Pk, 0)
}
if (logbase=="e"){ # same code as above, just with "e"
H_Nk <- apply(preyprop.mat, 2, function(x) -sum(x*log(x), na.rm=TRUE))
H_Pk <- apply(predprop.mat, 1, function(x) -sum(x*log(x), na.rm=TRUE))
n_Nk <- ifelse(colSums(web)!= 0, exp(H_Nk), 0)
n_Pk <- ifelse(rowSums(web)!= 0, exp(H_Pk), 0)
}
# mean number of predators per prey
V <- sum(rowSums(web)/sum(web)*n_Pk)
# marginal totals-weighted mean exp(Shannon diversity)
G <- sum(colSums(web)/sum(web)*n_Nk)
# linkage density
LD_q <- 0.5*(V+G)
#------------------
if ("linkage density" %in% index) out$"linkage density" <- LD_q
if ("weighted connectance" %in% index) out$"weighted connectance" <- LD_q/sum(dim(web))
#if ("weighted connectance" %in% index) out$"weighted connectance2" <- LD_q/mean(dim(web))
#LD_qs <- LD_q/(NROW(web)+NCOL(web)) # "weighted food web connectance", according to Jason's appendix
#------------------
# #We found no reference to this metric and saw little use for it. It is very similar to vulnerability/generality and can easily be computed from the output of \code{\link{specieslevel}} as \code{mean(specieslevel(web, index="diversity"))}
# if ("mean interaction diversity" %in% index){
# out$"interaction diversity LTL" <- mean(H_Nk)
# out$"interaction diversity HTL" <- mean(H_Pk)
#}
}
#------------------
if ("Fisher alpha" %in% index) {
fish <- try(fisherfit(web)$estimate, silent=TRUE) #in vegan
if (inherits(fish, "try-error")) {
out$"Fisher alpha" <- NA
} else {
out$"Fisher alpha" <- fish
}
}
#----------------------------- evenness & diversity ----------------------
if (any(c("interaction evenness", "Alatalo interaction evenness", "Shannon diversity") %in% index)){
# interaction evenness
p_i.mat <- web/sum(web)
#---------------
SH <- -sum(p_i.mat*log(p_i.mat), na.rm=TRUE)
if ("Shannon diversity" %in% index) out$"Shannon diversity" <- SH
IE <- ifelse(intereven=="prod", SH/log(prod(dim(web))), SH/log(sum(web>0)))
#---------------
if ("interaction evenness" %in% index) out$"interaction evenness" <- IE
#---------------
if ("Alatalo interaction evenness" %in% index){
evenness <- function(web){
# calculates evenness of the numbers of individual of different species in
# a community, NOT according to formula in Mueller et al. (1999,
# J. Anim. Ecol), but according to the original formula in Alatalo
# (1981, Oikos)
# can be extended at some point to more indices ...
pk <- web/sum(web)
(Alatalo <- (1/sum(pk^2) -1) / (exp(-sum(pk * log(pk), na.rm=TRUE)) -1))
}
E <- evenness(web)
out$"Alatalo interaction evenness" <- E
}
}
#---------------
# Bluethgen's H2'
if ("H2" %in% index){
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol # from the help of is.integer!
if (any(is.wholenumber(web)==FALSE)) H2_integer <- FALSE # turns H2_integer off if values are not integers
H2 <- as.numeric(H2fun(web, H2_integer=H2_integer)[1]) #1.element is the standardised H2 prime
out$"H2"= ifelse(H2<0, 0, H2)
}
#----------------------- now: grouplevel -------------------
# a list of network indices (which should not be called through grouplevel):
netw.index <- match(c("connectance", "web asymmetry", "links per species", "number of compartments", "compartment diversity", "modularity", "nestedness", "NODF", "weighted nestedness", "weighted NODF", "ISA", "SA", "interaction evenness", "Alatalo interaction evenness", "Fisher alpha", "H2", "Shannon diversity", "linkage density", "weighted connectance"), index)
exclude.index <- netw.index[!is.na(netw.index)]
gindex <- if (length(exclude.index)==0) index else index[-exclude.index] # exclude NAs from this vector
if (length(gindex) > 0) outg <- grouplevel(web, index=gindex, level=level, weighted=weighted, extinctmethod=extinctmethod, nrep=nrep, CCfun=CCfun, dist=dist, normalise=normalise, empty.web=empty.web, logbase=logbase, fcweighted=fcweighted, fcdist=fcdist)
if (exists("outg")){
# coerce potential list of HL/LL to one list:
if (is.list(outg)){
SEQ <- seq(1, 2*length(outg[[1]]), by=2)
sorted.outg <- c(outg[[1]], outg[[2]])
outg <- sorted.outg[order(c(SEQ, SEQ+1))]
}
out <- c(unlist(out), outg)
} else {out <- unlist(out)}
} # end legacy = FALSE
if (legacy == TRUE){ # for the old way of running networklevel
out <- .networklevel(web, index=index, ISAmethod=ISAmethod, SAmethod=SAmethod, extinctmethod=extinctmethod, nrep=nrep, plot.it.extinction=FALSE, plot.it.dd=FALSE, CCfun=CCfun, dist=dist, normalise=normalise, empty.web=empty.web, logbase=logbase, fcweighted=fcweighted, fcdist=fcdist)
} # end legacy condition
return(out)
}
#! JFedit: generality / vs. effective partners etc. NOT changed here! should it be? the changes are not really "legacy"?
.networklevel <- function(web, index="ALLBUTDD", ISAmethod="Bluethgen", SAmethod="Bluethgen", extinctmethod="r", nrep=100, plot.it.extinction=FALSE, plot.it.dd=FALSE, CCfun=median, dist="horn", normalise=TRUE, empty.web=TRUE, logbase="e", intereven="sum", H2_integer=TRUE, fcweighted=TRUE, fcdist="euclidean"){
##
## web interaction matrix, with lower trophic level in rows, higher in columns
## legacy se to TRUE allows to run networklevel in its old form
if(empty.web) {web <- empty(web)}
web.e <- empty(web) # emptied web for some indices
if (NROW(web) < 2 | NCOL(web) <2) warning("Web is really too small to calculate any reasonable index. You will get the values nonetheless, but I wouldn't put any faith in them!")
allindex <- c( #descriptive:
"number of species", "connectance", "web asymmetry",
#binary based and related:
"links per species", "number of compartments", "compartment diversity",
"cluster coefficient", "degree distribution", "mean number of shared partners", "modularity",
"togetherness", "C score", "V ratio", "discrepancy", "nestedness",
"weighted nestedness",
#miscelleneous:
"ISA", "SA", "extinction slope", "robustness", "niche overlap",
#quantitative series:
"weighted cluster coefficient", "weighted NODF", "generality", "vulnerability", "linkage density", "Fisher alpha", "interaction evenness", "Alatalo interaction evenness", "diversity", "effective partners", "functional complementarity", "H2" )
# GONE: "mean interaction diversity",
# only if indices are not given explicitly:
if (length(index)==1 & !all(index %in% allindex)){
index <- switch(index,
"ALL" = allindex,
"ALLBUTDD" = allindex[-which(allindex=="degree distribution")],
"info" = c("number of species", "connectance", "web asymmetry", "links per species", "number of compartments"),
# logic: only rough information on the network's general structure
"quantitative" = c("weighted cluster coefficient", "weighted nestedness", "weighted NODF", "functional complementarity", "H2", "diversity", "effective partners", "mean interaction diversity", "linkage density"),
# logic: the "quantitative series"
"binary" = c("connectance", "links per species", "nestedness", "cluster coefficient", "C-score"),
# logic: metrics for binary networks
"topology" = c("connectance", "cluster coefficient", "degree distribution", "togetherness", "nestedness"),
# logic: more abstract, topological metrics for binary networks
stop("Your index is not recognised! Typo? Check help for options!", call.=FALSE) #default for all non-matches
)
}
out <- list()
# set up enough panels for plotting:
if (plot.it.extinction) {m=2; n=1; par(mfrow=c(m,n), mar=c(5,5,4,1))} else m=1
if (plot.it.dd) {n=2; par(mfrow=c(m,n), mar=c(5,5,4,1))} else n=1
#-------------------
if ("number of species" %in% index) {
out$"number of lower trophic species" <- as.integer(NROW(web))
out$"number of higher trophic species" <- as.integer(NCOL(web))
}
#--------------------
if ("connectance" %in% index){
# connectance: "the fraction of all possible links that are realized in a network", (p. 12917 in Dunne et al. 2002)
out$connectance <- sum(web>0)/prod(dim(web))
}
#--------------------
if ("web asymmetry" %in% index) out$"web asymmetry" <- (NCOL(web)-NROW(web))/sum(dim(web)) # web asymmetry (Bluethgen et al. 2007, Fig. S2)
###---###---###---###---###---###---###---###---###---###---###---###---###
#--------------------
if ("links per species" %in% index){
L <- sum(web>0)/sum(dim(web))
out$"links per species"=L
}
#--------------------
if (any(c("number of compartments", "compartment diversity") %in% index)){
CD <- function(co){
if (co$n.compart>1){
no <- NA
for (i in 1:co$n.compart){
comp <- which(abs(co$cweb)==i, arr.ind=TRUE) # who is in a compartment?
no[i] <- length(unique(comp[,1])) + length(unique(comp[,2])) # how many species
}
no <- no/sum(dim(web)) # standardise for number of species in the web
CD <- exp(-sum(no*log(no)))
} else {CD <- NA} #; warning("only one compartment")}
CD
}
comps <- try(compart(web.e), silent=TRUE)
if (inherits(comps, "try-error")) {
ncompart <- compdiv <- NA
} else {
ncompart <- comps$n.compart
compdiv <- CD(comps)
}
if ("number of compartments" %in% index) out$"number of compartments" <- as.integer(ncompart)
if ("compartment diversity" %in% index) out$"compartment diversity" <- compdiv
}
#------------------------
if ("cluster coefficient" %in% index){
cluster.coef <- function(web, full=FALSE, FUN=mean){
# calculate cluster coefficient
# web a bipartite web
# full logical; return cluster coefficients for each species in the network?
# FUN give a function to summarise individual cluster coefficients, defaults to 'mean'.
# The concept was introduced by Watts & Strogatz (1998) and is described under in Wikipedia under http://en.wikipedia.org/w/index.php?title=Clustering_coefficient
# Its main idea was to help identifying small-world networks, which should
# have high cluster coefficients but low path lengths.
# Ci of species i is simply the number of realised links devided by the number of possible links. At the network level, the CC for a network is the average CC of its members. This is a little bit fishy, since the CCs are log-normal
# distributed (in pollinations networks at least). Therefore with FUN also
# other summary measures can be employed.
# Because within 'bipartite' we look at 2-mode networks, mean C is:
# C_lowertrophiclevel = C_highertrophiclevel = C_entirenetwork.
#
# Literature: Watts DJ, Strogatz S (1998) Collective dynmics of 'small-world' networks. Nature 393:440-442
# notice that this CC differs from the binary version of Tore Opsahl's. This here is the same as the one-mode cluster coefficient, Tore's is the "corrected" two-mode!
web <- as.matrix(web)
Ci.high <- colSums((web>0))/nrow(web)
Ci.low <- rowSums((web>0))/ncol(web)
CC <- FUN(Ci.high)
if (full) out <- list("cluster coefficient"=CC, "CC values higher"=Ci.high,
"CC values lower"=Ci.low) else out <- c("cluster coefficient"=CC)
out
}
out$"cluster coefficient"=as.numeric(cluster.coef(web, FUN=CCfun, full=FALSE))
}
if ("weighted cluster coefficient" %in% index){
# compute the weighted cluster coefficient using tnet:
edgelist <- web2edges(web, return=TRUE)
wcc <- try(unname(clustering_tm(edgelist)["am"]), silent=TRUE) #uses arithmetic mean!
out$"weighted cluster coefficient" <- if (inherits(wcc, "try-error")) "NA" else wcc
}
#--------------------------
# degree distribution fits:
if ("degree distribution" %in% index){
dd <- suppressWarnings(try(degreedistr(web, plot.it=plot.it.dd, pure.call=FALSE), silent=TRUE))
if (inherits(dd, "try-error")){
dd <- list()
dd$"dd fits LTL" <- NA
dd$"dd fits HTL" <- NA
}
out$"degree distribution LTL" <- dd$"dd fits LTL"
out$"degree distribution HTL" <- dd$"dd fits HTL"
}
#-------------------
# The Stone & Roberts's indices: S, T and C score:
#-------------------
if ("mean number of shared partners" %in% index) {
out$"mean number of shared HTL species" <-
mean(designdist(t(web)>0, method="J", terms="minimum"))
out$"mean number of shared LTL species" <-
mean(designdist(web>0, method="J", terms="minimum"))
}
#-------------------
if ("togetherness" %in% index){
out$togetherness <- togetherness(web, normalise=normalise, na.rm=TRUE)
}
#-------------------
if ("C score" %in% index){
out$"C score" <- C.score(web, normalise=normalise, na.rm=TRUE)
}
#-------------------
if ("V ratio" %in% index){
out$"V ratio" <- V.ratio(web)
}
#-------------------
if ("discrepancy" %in% index){
out$discrepancy <- as.integer(unname(vegan::nesteddisc(web)$statistic))
}
#-------------------
if ("nestedness" %in% index){
nest <- try(nestedtemp(web)$statistic, silent=TRUE)
out$nestedness <- ifelse(inherits(nest, "try-error"), NA, nest)
# a fast implementation of nestedness by Jari Oksanen
#old: nestedness(web, null.models=FALSE)$temperature
}
#-------------------
if ("weighted nestedness" %in% index){
wine.res <- try(wine(web.e, nreps=nrep)$wine, silent=TRUE)
out$"weighted nestedness" <- if (!inherits(wine.res, "try-error")) {wine.res} else {NA}
}
#-------------------
if ("weighted NODF" %in% index){
out$"weighted NODF" <- unname(nestednodf(web, order=TRUE, weighted=TRUE)$statistic[3])
}
###---###---###---###---###---###---###---###---###---###---###---###---###
if ("ISA" %in% index){#(any(c("ISA", "interaction strength asymmetry", "dependence asymmetry")) %in% index){
# Dependence asymmetry (Bascompte et al. 2006; Bluethgen et al. 2007, Fig. S2)
depL <- web.e/matrix(rowSums(web.e), nrow=NROW(web.e), ncol=NCOL(web.e), byrow=FALSE)
depH <- web.e/matrix(colSums(web.e), nrow=NROW(web.e), ncol=NCOL(web.e), byrow=TRUE)
if (ISAmethod=="Bascompte" & "ISA" %in% index) {
#depMax <- depL
#greaterindex <- depL < depH
#depMax[greaterindex] <- depH[greaterindex]
out$"dependence asymmetry"=mean(abs(depL-depH)/pmax(depL, depH), na.rm=TRUE)
}
if (ISAmethod=="Bluethgen" & "ISA" %in% index) {
web2 <- web
# delete cells for species encountered only once (Bluethgen, pers. comm.):
web2[, which(colSums(web)==1)] <- 0
web2[which(rowSums(web)==1), ] <- 0
rowsummat <- matrix(rowSums(web2), nrow=NROW(web2), ncol=NCOL(web2), byrow=FALSE)
colsummat <- matrix(colSums(web2), nrow=NROW(web2), ncol=NCOL(web2), byrow=TRUE)
depL <- web2/rowsummat
depH <- web2/colsummat
depL[depL<=0] <- NA
depH[depH<=0] <- NA
# now we need a correction to account for the fact that links with few (e.g. 2) observations will have a minimum depL of 1/2: all on one species: depL=1, one on each of two: depL=0.5
depLprime <- (depL - 1/rowsummat)/(1 - 1/rowsummat)
# assumes depLmin = 1/web2 and depLmax=1
depHprime <- (depH - 1/colsummat)/(1 - 1/colsummat)
out$"interaction strength asymmetry"=mean(as.matrix(depHprime-depLprime), na.rm=TRUE) #ranges from -1 to 1 /sum(depLprime, depHprime, na.rm=TRUE)
}
}
#---------------------------------------------------------
# Specialisation asymmetry (Bluethgen et al. 2007, Fig. S2)
# 2 options for calculating the "mean" SA:
# either as Bluethgen et al: average weighted by number of interactions in the
# cell or as mean of logarithms (since the dependencies follow a lognormal
# distribution)
if ("SA" %in% index){#(any(c("SA", "specialisation asymmetry")) %in% index){
di <- dfun(web)$dprime # plants
dj <- dfun(t(web))$dprime # pollinators
if (SAmethod=="log"){
lgmeani <- mean(log(di[di>0])); lgmeanj <- mean(log(dj[dj>0]))
SA <- (lgmeanj-lgmeani)/sum(lgmeani, lgmeanj)
# ij-sequence changed because log changes sequence, too
}
if (SAmethod=="Bluethgen"){
wmeani <- sum(di*rowSums(web.e))/sum(web.e)
wmeanj <- sum(dj*colSums(web.e))/sum(web.e)
SA <- (wmeanj-wmeani)/sum(wmeani, wmeanj)
# positive values indicate more specialisation in the higher trophic level
}
out$"specialisation asymmetry" <- SA
}
#--------------------------
# species extinction curve:
if (any(c("extinction slope", "robustness") %in% index)){
extL <- try(second.extinct(web=web, method=extinctmethod, nrep=nrep, participant="lower"), silent=TRUE)
extH <- try(second.extinct(web=web, method=extinctmethod, nrep=nrep, participant="higher"), silent=TRUE)
if ("extinction slope" %in% index){
slopeL <- try(slope.bipartite(extL, col="green", pch=16, type="b", plot.it=plot.it.extinction), silent=TRUE)
out$"extinction slope LTL"=suppressWarnings(as.numeric(slopeL))
slopeH <- try(slope.bipartite(extH, col="green", pch=16, type="b", plot.it=plot.it.extinction), silent=TRUE)
out$"extinction slope HTL"=suppressWarnings(as.numeric(slopeH))
}
if ("robustness" %in% index) {
if (inherits(extL, "try-error")){
robustL <- robustH <- NA
} else {
robustL <- try(robustness(extL), silent=TRUE)
robustH <- try(robustness(extH), silent=TRUE)
rL <- if (inherits(robustL, "try-error")) NA else robustL
rH <- if (inherits(robustH, "try-error")) NA else robustH
}
out$"robustness higher exterminated" = as.numeric(rH)
out$"robustness lower exterminated" = as.numeric(rL)
}
}
#--------------------------
# mean similarity of niches (niche overlap, sensu Krebs, Ecological Methodology)
# vegdist demands "sites" to be in rows, therefore the web has to be transposed
# to calculate dissimilarity between higher level species; similarity is simply
# 1-dissimilarity:
if ("niche overlap" %in% index) {
NOhigher <- mean(1-vegdist(t(web.e), method=dist))
NOlower <- mean(1-vegdist(web.e, method=dist))
out$"niche overlap LTL" <- NOlower
out$"niche overlap HTL" <- NOhigher
}
###---###---###---###---###---###---###---###---###---###---###---###---###
if (any(c("links per species", "linkage density", "vulnerability", "generality", "Fisher alpha", "mean interaction diversity") %in% index)){
# for formula see Tylianakis et al. (2006), supplement.
# N refers to prey, P to predators
preytot.mat <- matrix(rep(colSums(web), NROW(web)), NROW(web), byrow=TRUE)
preyprop.mat <- web/preytot.mat # = b_ik/b_.k in the first formula
#H_Nk is the diversity index of inflow (diversity of flower visits for each pollinator)
predtot.mat <- matrix(rep(rowSums(web), NCOL(web)), NROW(web), byrow=FALSE)
predprop.mat <- web/predtot.mat # = b_kj/b_.k in the second formula
if (logbase==2 | logbase=="2"){
H_Nk <- apply(preyprop.mat, 2, function(x) -sum(x*log2(x), na.rm=TRUE))
#H_Pk is the diversity index of pollinators for each plant species
H_Pk <- apply(predprop.mat, 1, function(x) -sum(x*log2(x), na.rm=TRUE))
# next, we need the reciprocals of this
# note that the ifelse is only needed if the web contains prey that is
# not eaten or predators that don't eat ...
n_Nk <- ifelse(colSums(web)!=0, 2^H_Nk, 0)
n_Pk <- ifelse(rowSums(web)!=0, 2^H_Pk, 0)
}
if (logbase=="e"){ # same code as above, just with "e"
H_Nk <- apply(preyprop.mat, 2, function(x) -sum(x*log(x), na.rm=TRUE))
H_Pk <- apply(predprop.mat, 1, function(x) -sum(x*log(x), na.rm=TRUE))
n_Nk <- ifelse(colSums(web)!= 0, exp(H_Nk), 0)
n_Pk <- ifelse(rowSums(web)!= 0, exp(H_Pk), 0)
}
# mean number of predators per prey
V <- sum(rowSums(web)/sum(web)*n_Pk)
# marginal totals-weighted mean exp(Shannon diversity)
G <- sum(colSums(web)/sum(web)*n_Nk)
# linkage density
LD_q <- 0.5*(V+G)
if ("vulnerability" %in% index) out$"vulnerability" <- V
#------------------
if ("generality" %in% index) out$"generality"=G
#------------------
if ("linkage density" %in% index) out$"linkage density"=LD_q
#LD_qs <- LD_q/(NROW(web)+NCOL(web)) # "weighted food web connectance", according to Jason's appendix
#------------------
if ("Fisher alpha" %in% index) {
fish <- try(fisherfit(web)$estimate, silent=TRUE) #in vegan
if (inherits(fish, "try-error")) {
out$"Fisher alpha" <- NA
} else {
out$"Fisher alpha" <- fish
}
}
#------------------
# #We found no reference to this metric and saw little use for it. It is very similar to vulnerability/generality and can easily be computed from the output of \code{\link{specieslevel}} as \code{mean(specieslevel(web, index="diversity"))}
# if ("mean interaction diversity" %in% index){
# out$"interaction diversity LTL" <- mean(H_Nk)
# out$"interaction diversity HTL" <- mean(H_Pk)
#}
}
#----------------------------- evenness & diversity ----------------------
if (any(c("interaction evenness", "Alatalo interaction evenness", "Shannon diversity") %in% index)){
# interaction evenness
p_i.mat <- web/sum(web)
SH <- -sum(p_i.mat*log(p_i.mat), na.rm=TRUE)
IE <- ifelse(intereven=="prod", SH/log(prod(dim(web))), SH/log(sum(web>0)))
#---------------
if ("interaction evenness" %in% index) out$"interaction evenness" <- IE
#---------------
if ("Alatalo interaction evenness" %in% index){
evenness <- function(web){
# calculates evenness of the numbers of individual of different species in
# a community, NOT according to formula in Mueller et al. (1999,
# J. Anim. Ecol), but according to the original formula in Alatalo
# (1981, Oikos)
# can be extended at some point to more indices ...
pk <- web/sum(web)
(Alatalo <- (1/sum(pk^2) -1) / (exp(-sum(pk * log(pk), na.rm=TRUE)) -1))
}
E <- evenness(web)
out$"Alatalo interaction evenness" <- E
}
#---------------
if ("Shannon diversity" %in% index) out$"Shannon diversity" <- SH
}
#---------------
# Devoto's fd:
#! JFedit: changed fd to fc here (but only the function and arguments because in legacy function; this may not be the best decision)
if (any(c("fc", "functional complementarity") %in% index)){
out$"Functional complementarity LTL" <- fc(t(web), dist=fcdist, method="average", weighted=fcweighted)
out$"Functional complementarity HTL" <- fc(web, dist=fcdist, method="average", weighted=fcweighted)
}
#---------------
# Bluethgen's H2'
if ("H2" %in% index){
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol # from the help of is.integer!
if (any(is.wholenumber(web)==FALSE)) H2_integer <- FALSE # turns H2_integer off if values are not integers
H2 <- as.numeric(H2fun(web, H2_integer=H2_integer)[1]) #1.element is the standardised H2 prime
out$"H2"= ifelse(H2<0, 0, H2)
}
#---------------------------------------------------------------------------
if (!("degree distribution" %in% index)) out <- unlist(out)
return(out)
}
#networklevel(Safariland, index="H2", legacy=F)
#networklevel(Safariland, index="ALLBUTDD", legacy=F)
#networklevel(Safariland, index=c("Alatalo interaction evenness", "H2"), legacy=F)
#networklevel(Safariland, index=c("extinction slope", "Alatalo interaction evenness", "H2"), legacy=F)
#networklevel(Safariland, index=c("H2", "degree distribution" ,"ISA", "SA", "clustering coefficient", "partner diversity"))
#networklevel(Safariland, index=c("H2", "extinction slope" ,"ISA", "SA", "cluster coefficient", "partner diversity", "diversity"))
#networklevel(Safariland, index="ALL")
#networklevel(Safariland, index=c("ISA", "SA"), level="both", weighted=FALSE)
#networklevel(Safariland,level="higher")
Try the bipartite package in your browser
Any scripts or data that you put into this service are public.
bipartite documentation built on Oct. 19, 2022, 1:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .networklevel `networklevel`
`networklevel` <- function(web, index="ALLBUTDD", level="both", weighted=TRUE, ISAmethod="Bluethgen", SAmethod="Bluethgen", extinctmethod="r", nrep=100, CCfun=median, dist="horn", normalise=TRUE, empty.web=TRUE, logbase="e", intereven="prod", H2_integer=TRUE, fcweighted=TRUE, fcdist="euclidean", legacy=FALSE){
##
## web interaction matrix, with lower trophic level in rows, higher in columns
## legacy se to TRUE allows to run networklevel in its old form
#! JFedit: fddist and fdweighted replace throughout by fcdist and fcweighted
## !!!!! when adding a new index, make sure to also add it to the "exclude from grouplevel"-list !!!!!##
if(empty.web) {web <- empty(web)}
web.e <- empty(web) # emptied web for some indices
if (NROW(web) < 2 | NCOL(web) <2) warning("Web is really too small to calculate any reasonable index. You will get the values nonetheless, but I wouldn't put any faith in them!")
allindex <- c( #descriptive:
"number of species", "connectance", "web asymmetry",
#binary based and related:
"links per species", "number of compartments", "compartment diversity",
"cluster coefficient", "degree distribution", "mean number of shared partners",
"togetherness", "C score", "V ratio", "discrepancy", "nestedness", "NODF", "weighted nestedness",
#miscelleneous:
"ISA", "SA", "extinction slope", "robustness", "niche overlap",
#quantitative series:
"weighted cluster coefficient", "modularity","weighted NODF", "partner diversity", "generality", "vulnerability", "linkage density", "weighted connectance", "Fisher alpha", "interaction evenness", "Alatalo interaction evenness", "Shannon diversity", "functional complementarity", "H2" )
# GONE: "mean interaction diversity", "effective partners",
index <- unique(index) # enforces that each index name is used only once
wrong.name <- which(is.na(pmatch(index, c(allindex, "ALL", "ALLBUTDD", "info","quantitative", "binary", "topology", "networklevel"))))
if (length(wrong.name) > 0) stop("You selected an index that is not available: ", paste(index[wrong.name], collapse=", "))
# only if indices are not given explicitly:
if (length(index)==1 & !all(index %in% allindex)){
index <- switch(index,
"ALL" = allindex,
"ALLBUTDD" = allindex[-which(allindex=="degree distribution")],
"info" = c("number of species", "connectance", "web asymmetry", "links per species", "number of compartments"),
# logic: only rough information on the network's general structure
"quantitative" = c("weighted cluster coefficient", "modularity", "weighted nestedness", "weighted NODF", "functional complementarity", "partner diversity", "effective partners", "H2", "diversity","linkage density", "weighted connectance", "niche overlap"), #"mean interaction diversity",
# logic: the "quantitative series"
"binary" = c("connectance", "links per species", "nestedness", "mean number of partners","cluster coefficient", "C-score", "Fisher alpha"),
# logic: metrics for binary networks
"topology" = c("connectance", "cluster coefficient", "degree distribution", "togetherness", "nestedness", "NODF"),
# logic: more abstract, topological metrics for binary networks
"networklevel" = c("connectance", "web asymmetry", "links per species", "number of compartments", "compartment diversity", "cluster coefficient", "modularity", "nestedness", "NODF", "weighted NODF", "ISA", "SA", "linkage density", "Fisher alpha", "diversity", "interaction evenness", "Alatalo interaction evenness", "H2"),
# only the truly networky indices
stop("Your index is not recognised! Typo? Check help for options!", call.=FALSE) #default for all non-matches
)
}
if (legacy == FALSE){
out <- list()
# first: all the true networklevel indices
#--------------------
if ("connectance" %in% index){
# connectance: "the fraction of all possible links that are realized in a network", (p. 12917 in Dunne et al. 2002)
suppressWarnings(out$connectance <- sum(web>0)/prod(dim(web)))
}
#--------------------
if ("web asymmetry" %in% index) out$"web asymmetry" <- (NCOL(web)-NROW(web))/sum(dim(web)) # web asymmetry (Bluethgen et al. 2007, Fig. S2)
#--------------------
if ("links per species" %in% index){
L <- sum(web>0)/sum(dim(web))
out$"links per species"=L
}
#--------------------
if (any(c("number of compartments", "compartment diversity") %in% index)){
CD <- function(co){
if (co$n.compart>1){
no <- NA
for (i in 1:co$n.compart){
comp <- which(abs(co$cweb)==i, arr.ind=TRUE) # who is in a compartment?
no[i] <- length(unique(comp[,1])) + length(unique(comp[,2])) # how many species
}
no <- no/sum(dim(web)) # standardise for number of species in the web
CD <- exp(-sum(no*log(no)))
} else {CD <- NA} #; warning("only one compartment")}
CD
}
comps <- try(compart(web.e), silent=TRUE)
if (inherits(comps, "try-error")) {
ncompart <- compdiv <- NA
} else {
ncompart <- comps$n.compart
compdiv <- CD(comps)
}
if ("number of compartments" %in% index) out$"number of compartments" <- as.integer(ncompart)
if ("compartment diversity" %in% index) out$"compartment diversity" <- compdiv
}
#------------------------
if ("cluster coefficient" %in% index){
cluster.coef <- function(web, full=FALSE, FUN=mean){
# calculate cluster coefficient
# web a bipartite web
# full logical; return cluster coefficients for each species in the network?
# FUN give a function to summarise individual cluster coefficients, defaults to 'mean'.
# The concept was introduced by Watts & Strogatz (1998) and is described under in Wikipedia under http://en.wikipedia.org/w/index.php?title=Clustering_coefficient
# Its main idea was to help identifying small-world networks, which should
# have high cluster coefficients but low path lengths.
# Ci of species i is simply the number of realised links devided by the number of possible links. At the network level, the CC for a network is the average CC of its members. This is a little bit fishy, since the CCs are log-normal
# distributed (in pollinations networks at least). Therefore with FUN also
# other summary measures can be employed.
# Because within 'bipartite' we look at 2-mode networks, mean C is:
# C_lowertrophiclevel = C_highertrophiclevel = C_entirenetwork.
#
# Literature: Watts DJ, Strogatz S (1998) Collective dynmics of 'small-world' networks. Nature 393:440-442
# notice that this CC differs from the binary version of Tore Opsahl's. This here is the same as the one-mode cluster coefficient, Tore's is the "corrected" two-mode!
web <- as.matrix(web)
Ci.high <- colSums((web>0))/nrow(web)
Ci.low <- rowSums((web>0))/ncol(web)
CC <- FUN(Ci.high)
if (full) out <- list("cluster coefficient"=CC, "CC values higher"=Ci.high,
"CC values lower"=Ci.low) else out <- c("cluster coefficient"=CC)
out
}
out$"cluster coefficient"=as.numeric(cluster.coef(web, FUN=CCfun, full=FALSE))
}
#-------------------
if ("modularity" %in% index){
Q <- try(computeModules(web)@likelihood)
out$"modularity Q" <- ifelse(inherits(Q, "try-error"), NA, Q)
}
#-------------------
if ("nestedness" %in% index){
nest <- try(nestedtemp(web)$statistic, silent=TRUE)
out$nestedness <- ifelse(inherits(nest, "try-error"), NA, nest)
# a fast implementation of nestedness by Jari Oksanen
#old: nestedness(web, null.models=FALSE)$temperature
}
#-------------------
if ("NODF" %in% index){
NODF <- try(unname(nestednodf(web, order=TRUE, weighted=FALSE)$statistic[3]), silent=TRUE)
out$"NODF" <- if (inherits(NODF, "try-error")) NA else NODF
}
#-------------------
if ("weighted nestedness" %in% index){
wine.res <- try(wine(web.e, nreps=nrep)$wine, silent=TRUE)
out$"weighted nestedness" <- if (!inherits(wine.res, "try-error")) {wine.res} else {NA}
}
#-------------------
if ("weighted NODF" %in% index){
wNODF <- try(unname(nestednodf(web, order=TRUE, weighted=TRUE)$statistic[3]), silent=TRUE)
out$"weighted NODF" <- if (inherits(wNODF, "try-error")) NA else wNODF
}
#------------------
if (any(c("ISA", "interaction strength asymmetry", "dependence asymmetry") %in% index)){
# Dependence asymmetry (Bascompte et al. 2006; Bluethgen et al. 2007, Fig. S2)
depL <- web.e/matrix(rowSums(web.e), nrow=NROW(web.e), ncol=NCOL(web.e), byrow=FALSE)
depH <- web.e/matrix(colSums(web.e), nrow=NROW(web.e), ncol=NCOL(web.e), byrow=TRUE)
if (ISAmethod=="Bascompte" & "ISA" %in% index) {
#depMax <- depL
#greaterindex <- depL < depH
#depMax[greaterindex] <- depH[greaterindex]
out$"dependence asymmetry"=mean(abs(depL-depH)/pmax(depL, depH), na.rm=TRUE)
}
if (ISAmethod=="Bluethgen" & "ISA" %in% index) {
web2 <- web
# delete cells for species encountered only once (Bluethgen, pers. comm.):
web2[, which(colSums(web)==1)] <- 0
web2[which(rowSums(web)==1), ] <- 0
rowsummat <- matrix(rowSums(web2), nrow=NROW(web2), ncol=NCOL(web2), byrow=FALSE)
colsummat <- matrix(colSums(web2), nrow=NROW(web2), ncol=NCOL(web2), byrow=TRUE)
depL <- web2/rowsummat
depH <- web2/colsummat
depL[depL<=0] <- NA
depH[depH<=0] <- NA
# now we need a correction to account for the fact that links with few (e.g. 2) observations will have a minimum depL of 1/2: all on one species: depL=1, one on each of two: depL=0.5
depLprime <- (depL - 1/rowsummat)/(1 - 1/rowsummat)
# assumes depLmin = 1/web2 and depLmax=1
depHprime <- (depH - 1/colsummat)/(1 - 1/colsummat)
out$"interaction strength asymmetry"=mean(as.matrix(depHprime-depLprime), na.rm=TRUE) #ranges from -1 to 1 /sum(depLprime, depHprime, na.rm=TRUE)
}
}
#------------------
if ("SA" %in% index){
# Specialisation asymmetry (Bluethgen et al. 2007, Fig. S2)
# 2 options for calculating the "mean" SA:
# either as Bluethgen et al: average weighted by number of interactions in the
# cell or as mean of logarithms (since the dependencies follow a lognormal
# distribution)
di <- dfun(web)$dprime # plants
dj <- dfun(t(web))$dprime # pollinators
if (SAmethod=="log"){
lgmeani <- mean(log(di[di>0])); lgmeanj <- mean(log(dj[dj>0]))
SA <- (lgmeanj-lgmeani)/sum(lgmeani, lgmeanj)
# ij-sequence changed because log changes sequence, too
}
if (SAmethod=="Bluethgen"){
wmeani <- sum(di*rowSums(web.e))/sum(web.e)
wmeanj <- sum(dj*colSums(web.e))/sum(web.e)
SA <- (wmeanj-wmeani)/sum(wmeani, wmeanj)
# positive values indicate more specialisation in the higher trophic level
}
out$"specialisation asymmetry" <- SA
}
#------------------
if (any(c("linkage density","weighted connectance") %in% index)){
# for formula see Tylianakis et al. (2006), supplement.
# N refers to prey, P to predators
preytot.mat <- matrix(rep(colSums(web), NROW(web)), NROW(web), byrow=TRUE)
preyprop.mat <- web/preytot.mat # = b_ik/b_.k in the first formula
#H_Nk is the diversity index of inflow (diversity of flower visits for each pollinator)
predtot.mat <- matrix(rep(rowSums(web), NCOL(web)), NROW(web), byrow=FALSE)
predprop.mat <- web/predtot.mat # = b_kj/b_.k in the second formula
if (logbase==2 | logbase=="2"){
H_Nk <- apply(preyprop.mat, 2, function(x) -sum(x*log2(x), na.rm=TRUE))
#H_Pk is the diversity index of pollinators for each plant species
H_Pk <- apply(predprop.mat, 1, function(x) -sum(x*log2(x), na.rm=TRUE))
# next, we need the reciprocals of this
# note that the ifelse is only needed if the web contains prey that is
# not eaten or predators that don't eat ...
n_Nk <- ifelse(colSums(web)!=0, 2^H_Nk, 0)
n_Pk <- ifelse(rowSums(web)!=0, 2^H_Pk, 0)
}
if (logbase=="e"){ # same code as above, just with "e"
H_Nk <- apply(preyprop.mat, 2, function(x) -sum(x*log(x), na.rm=TRUE))
H_Pk <- apply(predprop.mat, 1, function(x) -sum(x*log(x), na.rm=TRUE))
n_Nk <- ifelse(colSums(web)!= 0, exp(H_Nk), 0)
n_Pk <- ifelse(rowSums(web)!= 0, exp(H_Pk), 0)
}
# mean number of predators per prey
V <- sum(rowSums(web)/sum(web)*n_Pk)
# marginal totals-weighted mean exp(Shannon diversity)
G <- sum(colSums(web)/sum(web)*n_Nk)
# linkage density
LD_q <- 0.5*(V+G)
#------------------
if ("linkage density" %in% index) out$"linkage density" <- LD_q
if ("weighted connectance" %in% index) out$"weighted connectance" <- LD_q/sum(dim(web))
#if ("weighted connectance" %in% index) out$"weighted connectance2" <- LD_q/mean(dim(web))
#LD_qs <- LD_q/(NROW(web)+NCOL(web)) # "weighted food web connectance", according to Jason's appendix
#------------------
# #We found no reference to this metric and saw little use for it. It is very similar to vulnerability/generality and can easily be computed from the output of \code{\link{specieslevel}} as \code{mean(specieslevel(web, index="diversity"))}
# if ("mean interaction diversity" %in% index){
# out$"interaction diversity LTL" <- mean(H_Nk)
# out$"interaction diversity HTL" <- mean(H_Pk)
#}
}
#------------------
if ("Fisher alpha" %in% index) {
fish <- try(fisherfit(web)$estimate, silent=TRUE) #in vegan
if (inherits(fish, "try-error")) {
out$"Fisher alpha" <- NA
} else {
out$"Fisher alpha" <- fish
}
}
#----------------------------- evenness & diversity ----------------------
if (any(c("interaction evenness", "Alatalo interaction evenness", "Shannon diversity") %in% index)){
# interaction evenness
p_i.mat <- web/sum(web)
#---------------
SH <- -sum(p_i.mat*log(p_i.mat), na.rm=TRUE)
if ("Shannon diversity" %in% index) out$"Shannon diversity" <- SH
IE <- ifelse(intereven=="prod", SH/log(prod(dim(web))), SH/log(sum(web>0)))
#---------------
if ("interaction evenness" %in% index) out$"interaction evenness" <- IE
#---------------
if ("Alatalo interaction evenness" %in% index){
evenness <- function(web){
# calculates evenness of the numbers of individual of different species in
# a community, NOT according to formula in Mueller et al. (1999,
# J. Anim. Ecol), but according to the original formula in Alatalo
# (1981, Oikos)
# can be extended at some point to more indices ...
pk <- web/sum(web)
(Alatalo <- (1/sum(pk^2) -1) / (exp(-sum(pk * log(pk), na.rm=TRUE)) -1))
}
E <- evenness(web)
out$"Alatalo interaction evenness" <- E
}
}
#---------------
# Bluethgen's H2'
if ("H2" %in% index){
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol # from the help of is.integer!
if (any(is.wholenumber(web)==FALSE)) H2_integer <- FALSE # turns H2_integer off if values are not integers
H2 <- as.numeric(H2fun(web, H2_integer=H2_integer)[1]) #1.element is the standardised H2 prime
out$"H2"= ifelse(H2<0, 0, H2)
}
#----------------------- now: grouplevel -------------------
# a list of network indices (which should not be called through grouplevel):
netw.index <- match(c("connectance", "web asymmetry", "links per species", "number of compartments", "compartment diversity", "modularity", "nestedness", "NODF", "weighted nestedness", "weighted NODF", "ISA", "SA", "interaction evenness", "Alatalo interaction evenness", "Fisher alpha", "H2", "Shannon diversity", "linkage density", "weighted connectance"), index)
exclude.index <- netw.index[!is.na(netw.index)]
gindex <- if (length(exclude.index)==0) index else index[-exclude.index] # exclude NAs from this vector
if (length(gindex) > 0) outg <- grouplevel(web, index=gindex, level=level, weighted=weighted, extinctmethod=extinctmethod, nrep=nrep, CCfun=CCfun, dist=dist, normalise=normalise, empty.web=empty.web, logbase=logbase, fcweighted=fcweighted, fcdist=fcdist)
if (exists("outg")){
# coerce potential list of HL/LL to one list:
if (is.list(outg)){
SEQ <- seq(1, 2*length(outg[[1]]), by=2)
sorted.outg <- c(outg[[1]], outg[[2]])
outg <- sorted.outg[order(c(SEQ, SEQ+1))]
}
out <- c(unlist(out), outg)
} else {out <- unlist(out)}
} # end legacy = FALSE
if (legacy == TRUE){ # for the old way of running networklevel
out <- .networklevel(web, index=index, ISAmethod=ISAmethod, SAmethod=SAmethod, extinctmethod=extinctmethod, nrep=nrep, plot.it.extinction=FALSE, plot.it.dd=FALSE, CCfun=CCfun, dist=dist, normalise=normalise, empty.web=empty.web, logbase=logbase, fcweighted=fcweighted, fcdist=fcdist)
} # end legacy condition
return(out)
}
#! JFedit: generality / vs. effective partners etc. NOT changed here! should it be? the changes are not really "legacy"?
.networklevel <- function(web, index="ALLBUTDD", ISAmethod="Bluethgen", SAmethod="Bluethgen", extinctmethod="r", nrep=100, plot.it.extinction=FALSE, plot.it.dd=FALSE, CCfun=median, dist="horn", normalise=TRUE, empty.web=TRUE, logbase="e", intereven="sum", H2_integer=TRUE, fcweighted=TRUE, fcdist="euclidean"){
##
## web interaction matrix, with lower trophic level in rows, higher in columns
## legacy se to TRUE allows to run networklevel in its old form
if(empty.web) {web <- empty(web)}
web.e <- empty(web) # emptied web for some indices
if (NROW(web) < 2 | NCOL(web) <2) warning("Web is really too small to calculate any reasonable index. You will get the values nonetheless, but I wouldn't put any faith in them!")
allindex <- c( #descriptive:
"number of species", "connectance", "web asymmetry",
#binary based and related:
"links per species", "number of compartments", "compartment diversity",
"cluster coefficient", "degree distribution", "mean number of shared partners", "modularity",
"togetherness", "C score", "V ratio", "discrepancy", "nestedness",
"weighted nestedness",
#miscelleneous:
"ISA", "SA", "extinction slope", "robustness", "niche overlap",
#quantitative series:
"weighted cluster coefficient", "weighted NODF", "generality", "vulnerability", "linkage density", "Fisher alpha", "interaction evenness", "Alatalo interaction evenness", "diversity", "effective partners", "functional complementarity", "H2" )
# GONE: "mean interaction diversity",
# only if indices are not given explicitly:
if (length(index)==1 & !all(index %in% allindex)){
index <- switch(index,
"ALL" = allindex,
"ALLBUTDD" = allindex[-which(allindex=="degree distribution")],
"info" = c("number of species", "connectance", "web asymmetry", "links per species", "number of compartments"),
# logic: only rough information on the network's general structure
"quantitative" = c("weighted cluster coefficient", "weighted nestedness", "weighted NODF", "functional complementarity", "H2", "diversity", "effective partners", "mean interaction diversity", "linkage density"),
# logic: the "quantitative series"
"binary" = c("connectance", "links per species", "nestedness", "cluster coefficient", "C-score"),
# logic: metrics for binary networks
"topology" = c("connectance", "cluster coefficient", "degree distribution", "togetherness", "nestedness"),
# logic: more abstract, topological metrics for binary networks
stop("Your index is not recognised! Typo? Check help for options!", call.=FALSE) #default for all non-matches
)
}
out <- list()
# set up enough panels for plotting:
if (plot.it.extinction) {m=2; n=1; par(mfrow=c(m,n), mar=c(5,5,4,1))} else m=1
if (plot.it.dd) {n=2; par(mfrow=c(m,n), mar=c(5,5,4,1))} else n=1
#-------------------
if ("number of species" %in% index) {
out$"number of lower trophic species" <- as.integer(NROW(web))
out$"number of higher trophic species" <- as.integer(NCOL(web))
}
#--------------------
if ("connectance" %in% index){
# connectance: "the fraction of all possible links that are realized in a network", (p. 12917 in Dunne et al. 2002)
out$connectance <- sum(web>0)/prod(dim(web))
}
#--------------------
if ("web asymmetry" %in% index) out$"web asymmetry" <- (NCOL(web)-NROW(web))/sum(dim(web)) # web asymmetry (Bluethgen et al. 2007, Fig. S2)
###---###---###---###---###---###---###---###---###---###---###---###---###
#--------------------
if ("links per species" %in% index){
L <- sum(web>0)/sum(dim(web))
out$"links per species"=L
}
#--------------------
if (any(c("number of compartments", "compartment diversity") %in% index)){
CD <- function(co){
if (co$n.compart>1){
no <- NA
for (i in 1:co$n.compart){
comp <- which(abs(co$cweb)==i, arr.ind=TRUE) # who is in a compartment?
no[i] <- length(unique(comp[,1])) + length(unique(comp[,2])) # how many species
}
no <- no/sum(dim(web)) # standardise for number of species in the web
CD <- exp(-sum(no*log(no)))
} else {CD <- NA} #; warning("only one compartment")}
CD
}
comps <- try(compart(web.e), silent=TRUE)
if (inherits(comps, "try-error")) {
ncompart <- compdiv <- NA
} else {
ncompart <- comps$n.compart
compdiv <- CD(comps)
}
if ("number of compartments" %in% index) out$"number of compartments" <- as.integer(ncompart)
if ("compartment diversity" %in% index) out$"compartment diversity" <- compdiv
}
#------------------------
if ("cluster coefficient" %in% index){
cluster.coef <- function(web, full=FALSE, FUN=mean){
# calculate cluster coefficient
# web a bipartite web
# full logical; return cluster coefficients for each species in the network?
# FUN give a function to summarise individual cluster coefficients, defaults to 'mean'.
# The concept was introduced by Watts & Strogatz (1998) and is described under in Wikipedia under http://en.wikipedia.org/w/index.php?title=Clustering_coefficient
# Its main idea was to help identifying small-world networks, which should
# have high cluster coefficients but low path lengths.
# Ci of species i is simply the number of realised links devided by the number of possible links. At the network level, the CC for a network is the average CC of its members. This is a little bit fishy, since the CCs are log-normal
# distributed (in pollinations networks at least). Therefore with FUN also
# other summary measures can be employed.
# Because within 'bipartite' we look at 2-mode networks, mean C is:
# C_lowertrophiclevel = C_highertrophiclevel = C_entirenetwork.
#
# Literature: Watts DJ, Strogatz S (1998) Collective dynmics of 'small-world' networks. Nature 393:440-442
# notice that this CC differs from the binary version of Tore Opsahl's. This here is the same as the one-mode cluster coefficient, Tore's is the "corrected" two-mode!
web <- as.matrix(web)
Ci.high <- colSums((web>0))/nrow(web)
Ci.low <- rowSums((web>0))/ncol(web)
CC <- FUN(Ci.high)
if (full) out <- list("cluster coefficient"=CC, "CC values higher"=Ci.high,
"CC values lower"=Ci.low) else out <- c("cluster coefficient"=CC)
out
}
out$"cluster coefficient"=as.numeric(cluster.coef(web, FUN=CCfun, full=FALSE))
}
if ("weighted cluster coefficient" %in% index){
# compute the weighted cluster coefficient using tnet:
edgelist <- web2edges(web, return=TRUE)
wcc <- try(unname(clustering_tm(edgelist)["am"]), silent=TRUE) #uses arithmetic mean!
out$"weighted cluster coefficient" <- if (inherits(wcc, "try-error")) "NA" else wcc
}
#--------------------------
# degree distribution fits:
if ("degree distribution" %in% index){
dd <- suppressWarnings(try(degreedistr(web, plot.it=plot.it.dd, pure.call=FALSE), silent=TRUE))
if (inherits(dd, "try-error")){
dd <- list()
dd$"dd fits LTL" <- NA
dd$"dd fits HTL" <- NA
}
out$"degree distribution LTL" <- dd$"dd fits LTL"
out$"degree distribution HTL" <- dd$"dd fits HTL"
}
#-------------------
# The Stone & Roberts's indices: S, T and C score:
#-------------------
if ("mean number of shared partners" %in% index) {
out$"mean number of shared HTL species" <-
mean(designdist(t(web)>0, method="J", terms="minimum"))
out$"mean number of shared LTL species" <-
mean(designdist(web>0, method="J", terms="minimum"))
}
#-------------------
if ("togetherness" %in% index){
out$togetherness <- togetherness(web, normalise=normalise, na.rm=TRUE)
}
#-------------------
if ("C score" %in% index){
out$"C score" <- C.score(web, normalise=normalise, na.rm=TRUE)
}
#-------------------
if ("V ratio" %in% index){
out$"V ratio" <- V.ratio(web)
}
#-------------------
if ("discrepancy" %in% index){
out$discrepancy <- as.integer(unname(vegan::nesteddisc(web)$statistic))
}
#-------------------
if ("nestedness" %in% index){
nest <- try(nestedtemp(web)$statistic, silent=TRUE)
out$nestedness <- ifelse(inherits(nest, "try-error"), NA, nest)
# a fast implementation of nestedness by Jari Oksanen
#old: nestedness(web, null.models=FALSE)$temperature
}
#-------------------
if ("weighted nestedness" %in% index){
wine.res <- try(wine(web.e, nreps=nrep)$wine, silent=TRUE)
out$"weighted nestedness" <- if (!inherits(wine.res, "try-error")) {wine.res} else {NA}
}
#-------------------
if ("weighted NODF" %in% index){
out$"weighted NODF" <- unname(nestednodf(web, order=TRUE, weighted=TRUE)$statistic[3])
}
###---###---###---###---###---###---###---###---###---###---###---###---###
if ("ISA" %in% index){#(any(c("ISA", "interaction strength asymmetry", "dependence asymmetry")) %in% index){
# Dependence asymmetry (Bascompte et al. 2006; Bluethgen et al. 2007, Fig. S2)
depL <- web.e/matrix(rowSums(web.e), nrow=NROW(web.e), ncol=NCOL(web.e), byrow=FALSE)
depH <- web.e/matrix(colSums(web.e), nrow=NROW(web.e), ncol=NCOL(web.e), byrow=TRUE)
if (ISAmethod=="Bascompte" & "ISA" %in% index) {
#depMax <- depL
#greaterindex <- depL < depH
#depMax[greaterindex] <- depH[greaterindex]
out$"dependence asymmetry"=mean(abs(depL-depH)/pmax(depL, depH), na.rm=TRUE)
}
if (ISAmethod=="Bluethgen" & "ISA" %in% index) {
web2 <- web
# delete cells for species encountered only once (Bluethgen, pers. comm.):
web2[, which(colSums(web)==1)] <- 0
web2[which(rowSums(web)==1), ] <- 0
rowsummat <- matrix(rowSums(web2), nrow=NROW(web2), ncol=NCOL(web2), byrow=FALSE)
colsummat <- matrix(colSums(web2), nrow=NROW(web2), ncol=NCOL(web2), byrow=TRUE)
depL <- web2/rowsummat
depH <- web2/colsummat
depL[depL<=0] <- NA
depH[depH<=0] <- NA
# now we need a correction to account for the fact that links with few (e.g. 2) observations will have a minimum depL of 1/2: all on one species: depL=1, one on each of two: depL=0.5
depLprime <- (depL - 1/rowsummat)/(1 - 1/rowsummat)
# assumes depLmin = 1/web2 and depLmax=1
depHprime <- (depH - 1/colsummat)/(1 - 1/colsummat)
out$"interaction strength asymmetry"=mean(as.matrix(depHprime-depLprime), na.rm=TRUE) #ranges from -1 to 1 /sum(depLprime, depHprime, na.rm=TRUE)
}
}
#---------------------------------------------------------
# Specialisation asymmetry (Bluethgen et al. 2007, Fig. S2)
# 2 options for calculating the "mean" SA:
# either as Bluethgen et al: average weighted by number of interactions in the
# cell or as mean of logarithms (since the dependencies follow a lognormal
# distribution)
if ("SA" %in% index){#(any(c("SA", "specialisation asymmetry")) %in% index){
di <- dfun(web)$dprime # plants
dj <- dfun(t(web))$dprime # pollinators
if (SAmethod=="log"){
lgmeani <- mean(log(di[di>0])); lgmeanj <- mean(log(dj[dj>0]))
SA <- (lgmeanj-lgmeani)/sum(lgmeani, lgmeanj)
# ij-sequence changed because log changes sequence, too
}
if (SAmethod=="Bluethgen"){
wmeani <- sum(di*rowSums(web.e))/sum(web.e)
wmeanj <- sum(dj*colSums(web.e))/sum(web.e)
SA <- (wmeanj-wmeani)/sum(wmeani, wmeanj)
# positive values indicate more specialisation in the higher trophic level
}
out$"specialisation asymmetry" <- SA
}
#--------------------------
# species extinction curve:
if (any(c("extinction slope", "robustness") %in% index)){
extL <- try(second.extinct(web=web, method=extinctmethod, nrep=nrep, participant="lower"), silent=TRUE)
extH <- try(second.extinct(web=web, method=extinctmethod, nrep=nrep, participant="higher"), silent=TRUE)
if ("extinction slope" %in% index){
slopeL <- try(slope.bipartite(extL, col="green", pch=16, type="b", plot.it=plot.it.extinction), silent=TRUE)
out$"extinction slope LTL"=suppressWarnings(as.numeric(slopeL))
slopeH <- try(slope.bipartite(extH, col="green", pch=16, type="b", plot.it=plot.it.extinction), silent=TRUE)
out$"extinction slope HTL"=suppressWarnings(as.numeric(slopeH))
}
if ("robustness" %in% index) {
if (inherits(extL, "try-error")){
robustL <- robustH <- NA
} else {
robustL <- try(robustness(extL), silent=TRUE)
robustH <- try(robustness(extH), silent=TRUE)
rL <- if (inherits(robustL, "try-error")) NA else robustL
rH <- if (inherits(robustH, "try-error")) NA else robustH
}
out$"robustness higher exterminated" = as.numeric(rH)
out$"robustness lower exterminated" = as.numeric(rL)
}
}
#--------------------------
# mean similarity of niches (niche overlap, sensu Krebs, Ecological Methodology)
# vegdist demands "sites" to be in rows, therefore the web has to be transposed
# to calculate dissimilarity between higher level species; similarity is simply
# 1-dissimilarity:
if ("niche overlap" %in% index) {
NOhigher <- mean(1-vegdist(t(web.e), method=dist))
NOlower <- mean(1-vegdist(web.e, method=dist))
out$"niche overlap LTL" <- NOlower
out$"niche overlap HTL" <- NOhigher
}
###---###---###---###---###---###---###---###---###---###---###---###---###
if (any(c("links per species", "linkage density", "vulnerability", "generality", "Fisher alpha", "mean interaction diversity") %in% index)){
# for formula see Tylianakis et al. (2006), supplement.
# N refers to prey, P to predators
preytot.mat <- matrix(rep(colSums(web), NROW(web)), NROW(web), byrow=TRUE)
preyprop.mat <- web/preytot.mat # = b_ik/b_.k in the first formula
#H_Nk is the diversity index of inflow (diversity of flower visits for each pollinator)
predtot.mat <- matrix(rep(rowSums(web), NCOL(web)), NROW(web), byrow=FALSE)
predprop.mat <- web/predtot.mat # = b_kj/b_.k in the second formula
if (logbase==2 | logbase=="2"){
H_Nk <- apply(preyprop.mat, 2, function(x) -sum(x*log2(x), na.rm=TRUE))
#H_Pk is the diversity index of pollinators for each plant species
H_Pk <- apply(predprop.mat, 1, function(x) -sum(x*log2(x), na.rm=TRUE))
# next, we need the reciprocals of this
# note that the ifelse is only needed if the web contains prey that is
# not eaten or predators that don't eat ...
n_Nk <- ifelse(colSums(web)!=0, 2^H_Nk, 0)
n_Pk <- ifelse(rowSums(web)!=0, 2^H_Pk, 0)
}
if (logbase=="e"){ # same code as above, just with "e"
H_Nk <- apply(preyprop.mat, 2, function(x) -sum(x*log(x), na.rm=TRUE))
H_Pk <- apply(predprop.mat, 1, function(x) -sum(x*log(x), na.rm=TRUE))
n_Nk <- ifelse(colSums(web)!= 0, exp(H_Nk), 0)
n_Pk <- ifelse(rowSums(web)!= 0, exp(H_Pk), 0)
}
# mean number of predators per prey
V <- sum(rowSums(web)/sum(web)*n_Pk)
# marginal totals-weighted mean exp(Shannon diversity)
G <- sum(colSums(web)/sum(web)*n_Nk)
# linkage density
LD_q <- 0.5*(V+G)
if ("vulnerability" %in% index) out$"vulnerability" <- V
#------------------
if ("generality" %in% index) out$"generality"=G
#------------------
if ("linkage density" %in% index) out$"linkage density"=LD_q
#LD_qs <- LD_q/(NROW(web)+NCOL(web)) # "weighted food web connectance", according to Jason's appendix
#------------------
if ("Fisher alpha" %in% index) {
fish <- try(fisherfit(web)$estimate, silent=TRUE) #in vegan
if (inherits(fish, "try-error")) {
out$"Fisher alpha" <- NA
} else {
out$"Fisher alpha" <- fish
}
}
#------------------
# #We found no reference to this metric and saw little use for it. It is very similar to vulnerability/generality and can easily be computed from the output of \code{\link{specieslevel}} as \code{mean(specieslevel(web, index="diversity"))}
# if ("mean interaction diversity" %in% index){
# out$"interaction diversity LTL" <- mean(H_Nk)
# out$"interaction diversity HTL" <- mean(H_Pk)
#}
}
#----------------------------- evenness & diversity ----------------------
if (any(c("interaction evenness", "Alatalo interaction evenness", "Shannon diversity") %in% index)){
# interaction evenness
p_i.mat <- web/sum(web)
SH <- -sum(p_i.mat*log(p_i.mat), na.rm=TRUE)
IE <- ifelse(intereven=="prod", SH/log(prod(dim(web))), SH/log(sum(web>0)))
#---------------
if ("interaction evenness" %in% index) out$"interaction evenness" <- IE
#---------------
if ("Alatalo interaction evenness" %in% index){
evenness <- function(web){
# calculates evenness of the numbers of individual of different species in
# a community, NOT according to formula in Mueller et al. (1999,
# J. Anim. Ecol), but according to the original formula in Alatalo
# (1981, Oikos)
# can be extended at some point to more indices ...
pk <- web/sum(web)
(Alatalo <- (1/sum(pk^2) -1) / (exp(-sum(pk * log(pk), na.rm=TRUE)) -1))
}
E <- evenness(web)
out$"Alatalo interaction evenness" <- E
}
#---------------
if ("Shannon diversity" %in% index) out$"Shannon diversity" <- SH
}
#---------------
# Devoto's fd:
#! JFedit: changed fd to fc here (but only the function and arguments because in legacy function; this may not be the best decision)
if (any(c("fc", "functional complementarity") %in% index)){
out$"Functional complementarity LTL" <- fc(t(web), dist=fcdist, method="average", weighted=fcweighted)
out$"Functional complementarity HTL" <- fc(web, dist=fcdist, method="average", weighted=fcweighted)
}
#---------------
# Bluethgen's H2'
if ("H2" %in% index){
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol # from the help of is.integer!
if (any(is.wholenumber(web)==FALSE)) H2_integer <- FALSE # turns H2_integer off if values are not integers
H2 <- as.numeric(H2fun(web, H2_integer=H2_integer)[1]) #1.element is the standardised H2 prime
out$"H2"= ifelse(H2<0, 0, H2)
}
#---------------------------------------------------------------------------
if (!("degree distribution" %in% index)) out <- unlist(out)
return(out)
}
#networklevel(Safariland, index="H2", legacy=F)
#networklevel(Safariland, index="ALLBUTDD", legacy=F)
#networklevel(Safariland, index=c("Alatalo interaction evenness", "H2"), legacy=F)
#networklevel(Safariland, index=c("extinction slope", "Alatalo interaction evenness", "H2"), legacy=F)
#networklevel(Safariland, index=c("H2", "degree distribution" ,"ISA", "SA", "clustering coefficient", "partner diversity"))
#networklevel(Safariland, index=c("H2", "extinction slope" ,"ISA", "SA", "cluster coefficient", "partner diversity", "diversity"))
#networklevel(Safariland, index="ALL")
#networklevel(Safariland, index=c("ISA", "SA"), level="both", weighted=FALSE)
#networklevel(Safariland,level="higher")
Try the bipartite package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.