# Community functions that operate on trophic links
.ResolveTrophicLinksToNodeNames <- function(community, links)
{
# Returns a matrix with the columns 'resource' and 'consumer', both
# containing node names.
# links should be a matrix or data.frame with columns 'resource' and
# 'consumer'. If either column in links contains characters, they should be
# node names. If not characters, they should be integer node indices.
# The function does not check to make sure that the links are present in
# TLPS(community)
if(!is.Community(community)) stop('Not a Community')
# Get appropriate columns from links matrix
if(!is.null(colnames(links)))
{
stopifnot(all(is.element(c('resource','consumer'), colnames(links))))
links <- links[,c('resource','consumer'),drop=FALSE]
}
else
{
# Assume first two columns
stopifnot(ncol(links)>=2)
links <- links[,c(1,2),drop=FALSE]
}
node <- unname(NP(community, 'node'))
# Get resource and consumer as a node name
f <- function(values)
{
if(is.factor(values))
{
return (as.character(values))
}
else if(!is.character(values))
{
stopifnot(all(values>0 & values<=length(node)))
return (node[values])
}
else
{
return (values)
}
}
links <- cbind(resource=f(links[,1]), consumer=f(links[,2]))
rownames(links) <- NULL
bad <- !links %in% node
if(any(bad))
{
stop(paste('The names [', paste(unique(links[bad]), collapse=','),
'] are not nodes', sep=''))
}
return (links)
}
.ResolveTrophicLinksToRowIndices <- function(community, links)
{
# Returns a vector of integer indices in TLPS(community) for each row in
# links; NA for links that do not exist.
# links should be a matrix or data.frame with columns 'resource' and
# 'consumer'. If either column in links contains characters, they should be
# node names. If not characters, they should be integer node indices.
if(!is.Community(community)) stop('Not a Community')
links <- .ResolveTrophicLinksToNodeNames(community, links)
tlp <- TLPS(community)
match <- apply(links, 1, function (row)
{
l <- which(row[1]==tlp[,1] & row[2]==tlp[,2])
if(1==length(l))
{
return (l)
}
else if(0==length(l))
{
return (NA)
}
else
{
# If more than one row matches this link, something has gone really
# wrong.
stop(paste('Link from [', row[1], '] to [', row[2], '] appears ',
length(l), ' times in trophic links', sep=''))
}
})
return (match)
}
PredationMatrix <- function(community, weight=NULL)
{
# Returns a predation matrix. Columns and rows are named by node.
# Columns are consumers, rows are resource. If weight is NULL, elements
# will be either 0 (no trophic link) or 1 (trophic link).
# If not NULL, weight should be the names of a link property; non-zero
# elements of the matrix will contain the value of the weight property.
if(!is.Community(community)) stop('Not a Community')
tlp <- TLPS(community, link.properties=weight)
np <- NPS(community)
n.nodes <- nrow(np)
# Create predation matrix
# Casts to as.integer keep pm as a matrix of integer rather than
# numeric, which I think what we want.
pm <- matrix(as.integer(0), ncol=n.nodes, nrow=n.nodes)
colnames(pm) <- rownames(pm) <- np$node
if(!is.null(tlp))
{
stopifnot(all(tlp$resource %in% np$node))
stopifnot(all(tlp$consumer %in% np$node))
# Get numerical indices of resources and consumers
res <- sapply(tlp$resource, function(n) return (which(n==np$node)))
con <- sapply(tlp$consumer, function(n) return (which(n==np$node)))
pm[res + n.nodes * (con-1)] <- as.integer(1)
if(!is.null(weight))
{
pm[res + n.nodes * (con-1)] <- tlp[,weight]
}
}
return (pm)
}
InDegree <- TrophicGenerality <- NumberOfResources <- function(community)
{
# The number of resources of each node
tlps <- TLPS(community)
if(is.null(tlps))
{
res <- rep(0, NumberOfNodes(community))
names(res) <- unname(NP(community, 'node'))
return (res)
}
else
{
return (with(tlps, sapply(NP(community, 'node'),
function(n) sum(consumer==n))))
}
}
NormalisedTrophicGenerality <- function(community)
{
# The number of resources of each node, normalised with
# L/S. Williams and Martinez 2000 Nature.
return ( (1/LinkageDensity(community)) * TrophicGenerality(community))
}
OutDegree <- TrophicVulnerability <- NumberOfConsumers <- function(community)
{
# The number of consumers of each node
tlps <- TLPS(community)
if(is.null(tlps))
{
res <- rep(0, NumberOfNodes(community))
names(res) <- unname(NP(community, 'node'))
return (res)
}
else
{
return (with(tlps, sapply(NP(community, 'node'),
function(n) sum(resource==n))))
}
}
NormalisedTrophicVulnerability <- function(community)
{
# The number of consumers of each node, normalised with
# L/S. Williams and Martinez 2000 Nature.
return ( (1/LinkageDensity(community)) * TrophicVulnerability(community))
}
IsBasalNode <- function(community)
{
# TRUE for those nodes that consume no others, excluding
# themselves; will therefore return TRUE for nodes that are cannibalistic,
# have no resources and some consumers. Biologically, this should not
# happen, but mathematically, it should be accounted for.
community <- RemoveCannibalisticLinks(community)
return (0==NumberOfResources(community) & 0<NumberOfConsumers(community))
}
BasalNodes <- function(community)
{
# Names of nodes that consume no others
return (names(which(IsBasalNode(community))))
}
IsNonBasalNode <- function(community)
{
# TRUE for those nodes that are not basal
return (!IsBasalNode(community))
}
NonBasalNodes <- function(community)
{
# Names of nodes that consume at least one other
return (names(which(IsNonBasalNode(community))))
}
IsTopLevelNode <- function(community)
{
# TRUE for those node that are consumed by no others, excluding
# themselves; will therefore return TRUE for cannibalistic top level
# predators.
community <- RemoveCannibalisticLinks(community)
return (0==NumberOfConsumers(community) & 0<NumberOfResources(community))
}
TopLevelNodes <- function(community)
{
# Names of nodes that are consumed by no others, excluding themselves.
return (names(which(IsTopLevelNode(community))))
}
IsNonTopLevelNode <- function(community)
{
# Cannibalism ignored
return (!IsTopLevelNode(community))
}
NonTopLevelNodes <- function(community)
{
# Cannibalism ignored
return (names(which(IsNonTopLevelNode(community))))
}
IsIntermediateNode <- function(community)
{
# Neither basal nor top level. Cannibalism ignored.
community <- RemoveCannibalisticLinks(community)
return (0<NumberOfResources(community) & 0<NumberOfConsumers(community))
}
IntermediateNodes <- function(community)
{
# Neither basal nor top level
return (names(which(IsIntermediateNode(community))))
}
IsIsolatedNode <- function(community)
{
# TRUE for those node with no resource and no consumers
# Cannibalism ignored
community <- RemoveCannibalisticLinks(community)
return (0==NumberOfResources(community) & 0==NumberOfConsumers(community))
}
IsolatedNodes <- function(community)
{
# Nodes with no resources and no consumers
return (names(which(IsIsolatedNode(community))))
}
IsConnectedNode <- function(community)
{
# Cannibalism ignored
community <- RemoveCannibalisticLinks(community)
return (!IsIsolatedNode(community))
}
ConnectedNodes <- function(community)
{
return (names(which(IsConnectedNode(community))))
}
ResourcesByNode <- function(community)
{
# A named list of length NumberOfNodes(). Elements are vectors containing
# those nodes that are resources.
tlps <- TLPS(community)
if(is.null(tlps))
{
return (lapply(NP(community, 'node'),
function(n) vector(mode='character')))
}
else
{
return (with(tlps, lapply(NP(community, 'node'),
function(n) resource[consumer==n])))
}
}
ResourcesOfNodes <- function(community, nodes)
{
# Node should be either a vector of names or a vector of integer indices
if(!is.Community(community)) stop('Not a Community')
nodes <- .ResolveToNodeIndices(community, nodes)
r <- ResourcesByNode(community)
if(1==length(nodes))
{
# A vector of resources
return (r[[nodes]])
}
else
{
# A list of resource by node
return (r[nodes])
}
}
ConsumersByNode <- function(community)
{
# A named list of length NumberOfNodes(). Elements are vectors containing
# those nodes that are consumers.
if(!is.Community(community)) stop('Not a Community')
tlps <- TLPS(community)
if(is.null(tlps))
{
return (lapply(NP(community, 'node'),
function(n) vector(mode='character')))
}
else
{
return (with(tlps, lapply(NP(community, 'node'),
function(n) consumer[resource==n])))
}
}
ConsumersOfNodes <- function(community, nodes)
{
if(!is.Community(community)) stop('Not a Community')
nodes <- .ResolveToNodeIndices(community, nodes)
c <- ConsumersByNode(community)
if(1==length(nodes))
{
# A vector of consumers
return (c[[nodes]])
}
else
{
# A list of consumers by node
return (c[nodes])
}
}
ResourcesAndConsumersByNode <- function(community)
{
# A named list of length NumberOfNodes(). Elements are vectors containing
# those nodes that are either resources or consumers.
rc <- mapply(c, ResourcesByNode(community), ConsumersByNode(community),
SIMPLIFY=FALSE)
return (lapply(rc, unique))
}
NumberOfTrophicLinks <- function(community)
{
# The total number of links in the web.
tlps <- TLPS(community)
if(is.null(tlps)) return (0)
else return (nrow(tlps))
}
LinkageDensity <- function(community)
{
# The number of links per node
return (NumberOfTrophicLinks(community) / NumberOfNodes(community))
}
Degree <- function(community)
{
# A vector of length NumberOfNodes() containing the number of links
# for each node. Food web links are directed so cannibalistic links
# count twice.
return (InDegree(community) + OutDegree(community))
}
DegreeDistribution <- function(community, cumulative=FALSE)
{
# A vector of proportions of nodes with 0:max(degree) trophic
# links.
degree <- Degree(community)
dd <- tabulate(factor(degree, levels=0:max(degree)))
names(dd) <- 0:max(degree)
stopifnot(sum(dd)==NumberOfNodes(community))
dd <- dd / sum(dd)
if(cumulative)
{
dd <- rev(cumsum(rev(dd)))
}
return (dd)
}
DirectedConnectance <- function(community)
{
# Martinez 1992 Ecological Monographs
return (NumberOfTrophicLinks(community) / NumberOfNodes(community)^2)
}
IsCannibal <- function(community)
{
# TRUE if the node is a consumer of itself
tlps <- TLPS(community)
if(is.null(tlps))
{
res <- rep(FALSE, NumberOfNodes(community))
names(res) <- unname(NP(community, 'node'))
return (res)
}
else
{
return (with(tlps, sapply(NP(community, 'node'),
function(n) 1==sum(resource==n & consumer==n))))
}
}
Cannibals <- function(community)
{
# The names of those nodes which consume themselves
return (names(which(IsCannibal(community))))
}
FractionCannibalistic <- function(community)
{
n <- IsCannibal(community)
return (sum(n)/length(n))
}
.MatrixInversionTL <- function(community, weight.by, include.isolated)
{
# Returns a vector of trophic levels. If the trophic levels can not be
# computed then there is a problem with the community's network topology
# (see comment from Rich below), in which case all elements of the
# returned vector will be NA.
# Use the method described in Williams and Martinez (2004) Am Nat, Limits
# to Trophic Levels and Omnivory in Complex Food Webs: Theory and Data.
# Algorithm proposed by Levine 1980 J Theor Biol, Several measures of
# trophic structure applicable to complex food webs
# weight.by should be either NULL or the name of a trophic link property
# or function that gives diet fraction or similar.
# If weight.by is NULL, the returned values are 'prey-averaged TL'
# (Williams and Martinez (2004) Am Nat. If weight.by is not NULL,
# returned values are 'flow-based TL'.
# This method is quick and easy and accounts for flows through loops.
if(0==NumberOfTrophicLinks(community))
{
tl <- rep(NA, NumberOfNodes(community))
names(tl) <- unname(NP(community, 'node'))
}
else
{
# If requested by the caller, weighting is added to the matrix by
# PredationMatrix().
# In the absence of weighting / diet fraction information weight
# each link will be equally weighted.
pm <- PredationMatrix(community, weight=weight.by)
# The algorithm needs consumers on rows
W <- t(pm)
# Make each row sum to one.
rs <- rowSums(W)
W <- W / matrix(rs, ncol=ncol(W), nrow=nrow(W))
W[0==rs,] <- 0 # Fix NA resulting from div by zero
# Identity matrix
I <- diag(ncol(W))
# Invert the matrix. From Rich 2011-08-17: "If this matrix inversion
# fails, there is an important problem with the network topology.
# For a food web to be energetically feasible, every node must be
# connected to a basal node. When the inversion fails it is because
# there is at least one node that has no connection to a basal node."
result <- tryCatch(solve(I-W), error = function(e) e)
if('error' %in% class(result))
{
tl <- rep(NA, ncol(pm))
names(tl) <- colnames(pm)
}
else
{
# Returned vector is named by node
tl <- rowSums(result) # Isolated nodes have tl of 1
if(!include.isolated)
{
tl[IsolatedNodes(community)] <- NA
}
}
}
return (tl)
}
PreyAveragedTrophicLevel <- function(community, include.isolated=TRUE)
{
# Williams and Martinez (2004) Am Nat, p 460
# "Prey-averaged TL is equal to 1 + the mean TL of all the consumer's
# trophic resources TLj = 1 + sum{i=1 to s} l_{ij} TL_i / n_j
# where nj is the number of prey species in the diet of species j.
# This equation is equivalent to equation (1) [flow-based trophic level],
# with each nonzero link strength pij p 1/nj, which assumes that a consumer
# consumes all its prey species equally. "
return (.MatrixInversionTL(community, weight.by=NULL,
include.isolated=include.isolated))
}
FlowBasedTrophicLevel <- function(community, weight.by, include.isolated=TRUE)
{
# Williams and Martinez (2004) Am Nat, p 460
# Flow-based trophic level
stopifnot(!missing(weight.by))
return (.MatrixInversionTL(community, weight.by=weight.by,
include.isolated=include.isolated))
}
TrophicLevels <- function(community, weight.by=NULL, include.isolated=TRUE)
{
# Returns a matrix of either six or seven different measures of trophic
# level: the matrix-inversion method of PreyAveragedTrophicLevel,
# the matrix-inversion method of FlowBasedTrophicLevel (if weight is not
# NULL) and the five chain-position methods.
patl <- PreyAveragedTrophicLevel(community,
include.isolated=include.isolated)
fbtl <- NULL
if(!is.null(weight.by))
{
fbtl <- FlowBasedTrophicLevel(community, weight.by=weight.by,
include.isolated=include.isolated)
}
# The following measures of trophic level are computed using the
# position of nodes in each unique food chain.
chain.stats <- TrophicChainsStats(community)
S <- NumberOfNodes(community)
if(is.null(chain.stats))
{
# If no trophic links or if NumberOfTrophicLinks()>0 but chains is
# NULL (will happen if all links are cannibalistic), all trophic levels
# are NA
stl <- swtl <- ltl <- lwtl <- catl <- rep(NA, S)
}
else
{
# Williams and Martinez (2004) Am Nat, p 460
# "Shortest TL is equal to 1 + the shortest chain length from a
# consumer to a basal species."
stl <- apply(chain.stats$node.pos.counts, 1, function(row)
{
if(all(0==row)) NA
else min(which(0!=row))
})
# Williams and Martinez (2004) Am Nat, p 460
# "Short-weighted TL is the average of shortest TL and prey-averaged
# TL. This gives a measure biased toward shorter food chains."
swtl <- sapply(1:S, function(n) mean(c(stl[n],patl[n])))
# Williams and Martinez (2004) Am Nat, p 460
# "Longest TL is equal to 1 + the longest chain length from a consumer
# to a basal species."
ltl <- apply(chain.stats$node.pos.counts, 1, function(row)
{
if(all(0==row)) NA
else max(which(0!=row))
})
# Williams and Martinez (2004) Am Nat, p 460
# "Long-weighted TL is the average of longest TL and prey-averaged TL.
# This gives a measure biased toward longer food chains. "
lwtl <- sapply(1:S, function(n) mean(c(ltl[n], patl[n])))
# I think that Williams and Martinez's Chain-averaged TL is the same
# as Jonsson et al's Trophic Height:
# Williams and Martinez (2004) Am Nat, p 460
# "Chain-averaged TL is equal to 1 + the average chain length of all
# paths from a species to a basal species (Martinez 1991; Polis 1991;
# Fussman and Heber 2002)."
# Jonsson et al (2005) AER, p 8
# "The trophic position of a species in a food chain is 1 + the
# number of species preceding it in the ordered list of species in the
# chain... Trophic height is the average trophic position of a species
# in all food chains of which it is a part."
# A additional comment on flow-based vs chain-averaged TL
# Williams and Martinez (2004) Am Nat, p 460
# "Flow-based TL and prey-averaged TL are both computed using the
# matrix algebra method of Levine (1980) based on summing an infinite
# geometric series that includes the contributions from all loops.
# In contrast, the computation of chain-averaged TL maintains
# tractability in complex food webs by only passing through a loop
# once (Martinez 1991)."
catl <- apply(chain.stats$node.pos.counts, 1, function(row)
{
if(all(0==row)) NA
else weighted.mean(1:length(row), row)
})
if(include.isolated)
{
# Isolated species will have been assigned a TL of NA in the
# following measures. Set these to 1.
i <- IsIsolatedNode(community)
stl[i] <- swtl[i] <- ltl[i] <- lwtl[i] <- catl[i] <- 1
}
}
return (cbind(ShortestTL=stl,
ShortWeightedTL=swtl,
LongestTL=ltl,
LongWeightedTL=lwtl,
ChainAveragedTL=catl,
PreyAveragedTL=patl,
FlowBasedTL=fbtl))
}
# Convenience functions to access just one of the chain-averaged measures
ShortestTrophicLevel <- function(community, include.isolated=TRUE)
{
tl <- TrophicLevels(community, weight.by=NULL,
include.isolated=include.isolated)
res <- tl[,'ShortestTL']
names(res) <- rownames(tl)
return (res)
}
ShortWeightedTrophicLevel <- function(community, include.isolated=TRUE)
{
tl <- TrophicLevels(community, weight.by=NULL,
include.isolated=include.isolated)
res <- tl[,'ShortWeightedTL']
names(res) <- rownames(tl)
return (res)
}
LongestTrophicLevel <- function(community, include.isolated=TRUE)
{
tl <- TrophicLevels(community, weight.by=NULL,
include.isolated=include.isolated)
res <- tl[,'LongestTL']
names(res) <- rownames(tl)
return (res)
}
LongWeightedTrophicLevel <- function(community, include.isolated=TRUE)
{
tl <- TrophicLevels(community, weight.by=NULL,
include.isolated=include.isolated)
res <- tl[,'LongWeightedTL']
names(res) <- rownames(tl)
return (res)
}
ChainAveragedTrophicLevel <-
TrophicHeight <- function(community, include.isolated=TRUE)
{
tl <- TrophicLevels(community, weight.by=NULL,
include.isolated=include.isolated)
res <- tl[,'ChainAveragedTL']
names(res) <- rownames(tl)
return (res)
}
TrophicSpecies <- function(community, include.isolated=TRUE)
{
# Returns a vector containing the trophic species number of each
# node in the community. Nodes with identical sets of resources and
# consumers are assigned to the same trophic species.
# If include.isolated is TRUE, isolated node are assigned a trophic
# node. If include.isolated is FALSE isolated nodes have a trophic
# node of NA.
# library(cheddar)
# data(Benguela, BroadstoneStream, TL84, TL86, SkipwithPond,YthanEstuary)
# for(community in list(Benguela, BroadstoneStream, TL84, TL86,
# SkipwithPond, YthanEstuary))
# {
# print(community)
# print(system.time(ts <- TrophicSpecies(community)))
# }
# Benguela containing 29 nodes.
# user system elapsed
# 0.008 0.000 0.009
# Broadstone Stream containing 37 nodes.
# user system elapsed
# 0.008 0.000 0.008
# Tuesday Lake sampled in 1984 containing 56 nodes.
# user system elapsed
# 0.012 0.000 0.012
# Tuesday Lake sampled in 1986 containing 57 nodes.
# user system elapsed
# 0.012 0.000 0.012
# Skipwith Pond containing 37 nodes.
# user system elapsed
# 0.012 0.000 0.012
# Ythan Estuary containing 92 nodes.
# user system elapsed
# 0.032 0.000 0.030
if(!is.Community(community)) stop('Not a Community')
trophic.species = rep(NA, NumberOfNodes(community))
names(trophic.species) <- unname(NP(community, 'node'))
remaining <- 1:NumberOfNodes(community)
if(!include.isolated)
{
remaining <- setdiff(remaining, which(IsIsolatedNode(community)))
}
resources <- ResourcesByNode(community)
consumers <- ConsumersByNode(community)
trophic.species.number <- 1
while(length(remaining)>0)
{
# Consider the first node that has not yet been assigned a trophic
# species
node <- remaining[1]
# Remove from remaining and assign trophic species number
remaining <- setdiff(remaining, node)
trophic.species[node] <- trophic.species.number
for(t in remaining)
{
if(length(resources[[node]])==length(resources[[t]]) &&
length(consumers[[node]])==length(consumers[[t]]) &&
all(resources[[node]]==resources[[t]]) &&
all(consumers[[node]]==consumers[[t]]))
{
# Remove from remaining and assign trophic species number
remaining <- setdiff(remaining, t)
trophic.species[t] <- trophic.species.number
}
}
trophic.species.number <- trophic.species.number + 1
}
return (trophic.species)
}
TrophicLinksForNodes <- function(community, nodes, node.properties=NULL,
link.properties=NULL)
{
# Returns a matrix of links in-to and out-of the given nodes.
# nodes should be either node numbers or node names.
if(!is.Community(community)) stop('Not a Community')
stopifnot((is.character(nodes) && all(nodes %in% NP(community, 'node'))) ||
all(0<nodes & nodes<=NumberOfNodes(community)))
if(!is.character(nodes))
{
nodes <- unname(NP(community, 'node')[nodes])
}
tlp <- TLPS(community, node.properties=node.properties,
link.properties=link.properties)
links <- tlp[tlp$resource %in% nodes | tlp$consumer %in% nodes,,drop=FALSE]
rownames(links) <- NULL
return (links)
}
.DoLinksExist <- function(community, links)
{
# Returns TRUE or FALSE for each row in links.
# links should be either a matrix or a data.frame.
# links should contain column called resource and consumer.
# Values should be either node numbers or node names.
l <- .ResolveTrophicLinksToRowIndices(community, links)
return (!is.na(l))
}
RemoveCannibalisticLinks <- function(community,
title=paste(CP(community, 'title'),
'(cannibalistic links removed)'))
{
# Returns a copy of community with cannibalistic links removed
if(!is.Community(community)) stop('Not a Community')
tlp <- TLPS(community)
if(is.null(tlp))
{
return (community)
}
else
{
remove <- tlp[,1] == tlp[,2]
tlp <- tlp[!remove,,drop=FALSE]
properties <- CPS(community)
properties$title <- title
return (Community(NPS(community), trophic.links=tlp,
properties=properties))
}
}
FractionBasalNodes <- function(community)
{
n <- IsBasalNode(community)
return (length(which(n)) / length(n))
}
FractionNonBasalNodes <- function(community)
{
n <- IsNonBasalNode(community)
return (length(which(n)) / length(n))
}
FractionTopLevelNodes <- function(community)
{
n <- IsTopLevelNode(community)
return (length(which(n)) / length(n))
}
FractionNonTopLevelNodes <- function(community)
{
n <- IsNonTopLevelNode(community)
return (length(which(n)) / length(n))
}
FractionIntermediateNodes <- function(community)
{
n <- IsIntermediateNode(community)
return (length(which(n)) / length(n))
}
FractionIsolatedNodes <- function(community)
{
n <- IsIsolatedNode(community)
return (length(which(n)) / length(n))
}
FractionConnectedNodes <- function(community)
{
n <- IsConnectedNode(community)
return (length(which(n)) / length(n))
}
ShortestPaths <- function(community, weight.by=NULL)
{
# Returns a matrix of NumberOfNodes square.
# Dijkstra's algorithm is used to compute path lengths.
# library(cheddar)
# data(Benguela, BroadstoneStream, TL84, TL86, SkipwithPond,YthanEstuary)
# for(community in list(Benguela, BroadstoneStream, TL84, TL86,
# SkipwithPond, YthanEstuary))
# {
# print(community)
# print(system.time(y <- ShortestPaths(community)))
# }
# Benguela containing 29 nodes.
# user system elapsed
# 0.024 0.000 0.024
# Broadstone Stream containing 37 nodes.
# user system elapsed
# 0.020 0.000 0.018
# Tuesday Lake sampled in 1984 containing 56 nodes.
# user system elapsed
# 0.032 0.000 0.034
# Tuesday Lake sampled in 1986 containing 57 nodes.
# user system elapsed
# 0.032 0.000 0.032
# Skipwith Pond containing 37 nodes.
# user system elapsed
# 0.032 0.000 0.034
# Ythan Estuary containing 92 nodes.
# user system elapsed
# 0.060 0.000 0.067
if(!is.Community(community)) stop('Not a Community')
# Get the food web as two adjancency lists
consumers <- .CAdjacencyList(community, ConsumersByNode(community))
resources <- .CAdjacencyList(community, ResourcesByNode(community))
pm <- PredationMatrix(community, weight=weight.by)
n <- NumberOfNodes(community)
lengths <- as.double(rep(NA, n*n))
status <- as.integer(-1)
res <- .C('shortest_paths',
as.integer(consumers),
as.integer(length(consumers)),
as.integer(resources),
as.integer(length(resources)),
as.double(pm),
as.integer(n),
lengths=lengths,
status=status,
PACKAGE='cheddar',
NAOK=TRUE)
if(0==res$status)
{
# Take the subset of the matrix that we are interested in
sp <- res$lengths
dim(sp) <- c(n,n)
colnames(sp) <- rownames(sp) <- unname(NP(community, 'node'))
}
else if(1==res$status)
{
stop('Problem with an input parameter')
}
else
{
stop(paste('Unknown status [', res$status, ']', sep=''))
}
return (sp)
}
CharacteristicPathLength <- function(community)
{
# From Williams and Martinez 2002 PNAS:
# "...characteristic path length, D quantifies the average number of
# links necessary for information or an effect to propagate along
# the shortest paths between nodes in a network"
spl <- ShortestPaths(community)
spl[is.infinite(spl)] <- NA
return (mean(spl, na.rm=TRUE))
}
.SumPredationMatrixGap <- function(community, diet.gap)
{
# Stouffer et al (2006) PNAS. Delegate to the C implementation not for
# speed but for consistency.
if(!is.Community(community)) stop('Not a Community')
pm <- PredationMatrix(community)
if(!diet.gap)
{
pm <- t(pm)
}
n <- as.integer(ncol(pm))
sum_diet_gaps <- as.integer(0)
status <- as.integer(-1)
res <- .C('sum_diet_gaps',
as.integer(pm),
as.integer(n),
as.integer((1:n) - 1), # Row order, 0-indexed
sum_diet_gaps=sum_diet_gaps,
status=status)
if(-1==res$status)
{
stop('Unexpected error')
}
else if(0==res$status)
{
# C 0-index back to R 1-index
return (res$sum_diet_gaps)
}
else if(1==res$status)
{
stop('Problem with an input parameter')
}
else(1==res$status)
{
stop(paste('Unknown status [', res$status, ']', sep=''))
}
}
SumDietGaps <- function(community)
{
# Stouffer et al (2006) PNAS. Delegate to the C implementation not for
# speed but for consistency.
return (.SumPredationMatrixGap(community, diet.gap=TRUE))
}
SumConsumerGaps <- function(community)
{
# Zook et al (2011) J Theor Biol. Delegate to the C implementation not for
# speed but for consistency.
return (.SumPredationMatrixGap(community, diet.gap=FALSE))
}
.MinimisePredationMatrixGaps <- function(community, T.start, T.stop, c,
swaps.per.T, trace.anneal,
diet.gap, n, best)
{
# Using Stouffer et al's (2006) PNAS simulated annealing learning method.
# If diet.gap is TRUE then the sum diet gap is minimized otherwise the
# sum in each species' consumers is minimized (ZookEtAl2011JTheorEcol).
# library(cheddar)
# data(Benguela, BroadstoneStream, TL84, TL86, SkipwithPond,YthanEstuary)
# for(community in list(Benguela, BroadstoneStream, TL84, TL86,
# SkipwithPond,YthanEstuary))
# {
# print(community)
# print(system.time(i <- MinimiseSumDietGaps(community)))
# print(paste('Diet gap:',i$sum.gaps))
# }
# Benguela containing 29 nodes.
# user system elapsed
# 0.144 0.000 0.144
# [1] "Diet gap: 27"
# Broadstone Stream containing 37 nodes.
# user system elapsed
# 0.080 0.004 0.083
# [1] "Diet gap: 7"
# Tuesday Lake sampled in 1984 containing 56 nodes.
# user system elapsed
# 0.204 0.000 0.204
# [1] "Diet gap: 8"
# Tuesday Lake sampled in 1986 containing 57 nodes.
# user system elapsed
# 0.184 0.004 0.185
# [1] "Diet gap: 9"
# Skipwith Pond containing 37 nodes.
# user system elapsed
# 0.156 0.000 0.159
# [1] "Diet gap: 31"
# Ythan Estuary containing 92 nodes.
# user system elapsed
# 0.512 0.000 0.513
# [1] "Diet gap: 262"
if(!is.Community(community)) stop('Not a Community')
stopifnot(T.start>T.stop)
stopifnot(T.stop>0)
stopifnot(c>0 && c<1)
stopifnot(swaps.per.T>0)
if(is.null(TLPS(community)))
{
stop('The community has no trophic links')
}
stopifnot(n>0)
pm <- PredationMatrix(community)
if(!diet.gap)
{
pm <- t(pm)
}
F <- function()
{
# Out-params
sum.gaps <- as.integer(-1)
best <- as.integer(rep(NA, NumberOfNodes(community)))
status <- as.integer(-1)
res <- .C('minimise_sum_diet_gaps',
as.integer(pm),
as.integer(NumberOfNodes(community)),
as.double(T.start),
as.double(T.stop),
as.double(c),
as.integer(swaps.per.T),
as.integer(trace.anneal),
sum.gaps=sum.gaps,
best=best,
status=status,
PACKAGE='cheddar',
NAOK=TRUE)
if(-1==res$status)
{
stop('Unexpected error')
}
else if(0==res$status)
{
# C 0-index back to R 1-index
best <- res$best + 1
stopifnot(all(best>0 & best<=NumberOfNodes(community)))
reordered <- OrderCommunity(community, new.order=best,
title=paste(CP(community, 'title'),
'sum',
ifelse(diet.gap, 'diet', 'consumer'),
'gap minimised to',
res$sum.gaps))
return (list(sum.gaps=res$sum.gaps,
order=unname(NP(community, 'node'))[best],
reordered=reordered))
}
else if(1==res$status)
{
stop('Problem with an input parameter')
}
else(1==res$status)
{
stop(paste('Unknown status [', res$status, ']', sep=''))
}
}
res <- replicate(n, F(), simplify=FALSE)
res <- res[order(sapply(res, '[[', 'sum.gaps'))]
if(best || 1==n)
{
return (res[[1]])
}
else
{
return (res)
}
}
MinimiseSumDietGaps <- function(community, T.start=10, T.stop=0.1, c=0.9,
swaps.per.T=1000, trace.anneal=FALSE, n=1,
best=TRUE)
{
return (.MinimisePredationMatrixGaps(community, T.start, T.stop, c,
swaps.per.T, trace.anneal,
diet.gap=TRUE, n, best))
}
MinimiseSumConsumerGaps <- function(community, T.start=10, T.stop=0.1, c=0.9,
swaps.per.T=1000, trace.anneal=FALSE,
n=1, best=TRUE)
{
return (.MinimisePredationMatrixGaps(community, T.start, T.stop, c,
swaps.per.T, trace.anneal,
diet.gap=FALSE, n, best))
}
TrophicSimilarity <- function(community)
{
# Martinez 1991 Ecol. Monog. 61, 367-392 p 370.
# I = c/(a+b+c)
# where c = number of predators and prey common to the two taxa,
# a = number of predators and prey unique to one taxon, and b = number
# of predators and prey unique to one taxon. When the two taxa have the
# same set of predators and prey, I = 1, When the two taxa have no
# common predators or common prey, I = 0.
rc <- ResourcesAndConsumersByNode(community)
nodes <- names(rc)
I <- matrix(NA, ncol=length(nodes), nrow=length(nodes),
dimnames=list(nodes,nodes))
diag(I) <- 1
if(1<NumberOfNodes(community))
{
# Inner loop depends upon the community having more than one node
for(x in 1:(length(nodes)-1))
{
for(y in (x+1):length(nodes))
{
a <- length(setdiff(rc[[x]], rc[[y]]))
b <- length(setdiff(rc[[y]], rc[[x]]))
c <- length(intersect(rc[[x]], rc[[y]]))
I[x,y] <- I[y,x] <- c/(a+b+c)
}
}
}
I[IsolatedNodes(community),IsolatedNodes(community)] <- NA
stopifnot(all.equal(I, t(I)))
stopifnot(all(is.na(I)) || 1>=max(I, na.rm=TRUE))
stopifnot(1==diag(I)[-which(IsIsolatedNode(community))])
return (I)
}
MeanMaximumTrophicSimilarity <- function(community)
{
# Williams and Martinez 2000 Nature 404 180-182, p 181.
# MxSim = 1/S \sum_1^S max(s_{ij}) (i!=j)
I <- TrophicSimilarity(community)
diag(I) <- NA # Exclude diagonal
return (mean(apply(I, 2, max, na.rm=TRUE)))
}
IsOmnivore <- function(community, level=PreyAveragedTrophicLevel)
{
# Species that consume two or more species and have a non-integer
# trophic level
# Polis, G. A. Complex desert food webs: an empirical critique of food web
# theory. Am. Nat. 138, 123-155 (1991).
n.resources <- sapply(ResourcesByNode(community), length)
tl <- level(community)
return (n.resources>=2 & 0!=(tl %% 1))
}
Omnivores <- function(community, ...)
{
return (names(which(IsOmnivore(community, ...))))
}
FractionOmnivorous <- Omnivory <- function(community, ...)
{
# The fraction of species that are omnivores
return (length(Omnivores(community, ...))/NumberOfNodes(community))
}
NodeQuantitativeDescriptors <- function(community, weight)
{
# Bersier et al (2002) Ecology
if(!is.Community(community)) stop('Not a Community')
.RequireTrophicLinks(community)
# A common operation
vlog2v <- function(v) v*log2(v)
b <- PredationMatrix(community, weight)
bIn <- colSums(b)
bOut <- rowSums(b)
# Diversity of inflows and outflows
HN <- -rowSums(vlog2v(t(b)/bIn), na.rm=TRUE) # p 2397, eq 5
HP <- -rowSums(vlog2v(b/bOut), na.rm=TRUE) # p 2397, eq 6
# Equivalent numbers of prey and predators
nN <-2^HN # p 2397, eq 7
nN[0==bIn] <- 0
nP <-2^HP # p 2397, eq 8
nP[0==bOut] <- 0
d.prime <- nN / (nN + nP) # p 2397, eq 9
d <- bIn*nN / (bIn*nN + bOut*nP) # p 2397, eq 10
# p 2396:
# Trophic level of a taxon (Yodzis 1989); the characterization we adopt
# here is 'one plus the length of the longest chain from the focal taxon
# to a basal taxon.'
t <- LongestTrophicLevel(community)
res <- ResourcesByNode(community)
o.prime <- -1 + 2^sapply(1:NumberOfNodes(community), function(k)
{
# Compute the number of prey taxa located at each trophic level t
r <- res[[k]] # Resources of k
n <- table(t[r]) # The number of resources at each trophic level
return (-sum(vlog2v(n/sum(n))))
})
o <- -1 + 2^sapply(1:NumberOfNodes(community), function(k)
{
r <- res[[k]] # Resources of k
# Biomass going into k summed by trophic level of resource
bt <- tapply(b[r,k], t[r], sum)
return (-sum(vlog2v(bt/sum(bt))))
})
# Standardised quantitative G and V, p 2400, eq 28-29
g.prime <- nN*NumberOfNodes(community) / sum(nN)
g <- bIn*nN*NumberOfNodes(community) / sum(bIn*nN)
v.prime <- nP*NumberOfNodes(community) / sum(nP)
v <- bOut*nP*NumberOfNodes(community) / sum(bOut*nP)
# Bersier et al table 1
res <- cbind(NResources=NumberOfResources(community),
NConsumers=NumberOfConsumers(community),
bIn,
bOut,
nN,
nP,
d.prime,
d,
o.prime,
o,
g.prime,
g,
v.prime,
v)
rownames(res) <- unname(NP(community, 'node'))
return (res)
}
QuantitativeDescriptors <- function(community, weight, top.level.threshold=0.99)
{
# A common operation
vlog2v <- function(v) v*log2(v)
# Community-level descriptors are derived from the node-level decsriptors
# TODO prevent double-computation of chains
np <- NodeQuantitativeDescriptors(community, weight)
b <- PredationMatrix(community, weight)
sumb <- sum(b)
# Diversity of inflows and outflows
HN <- -rowSums(vlog2v(t(b)/np[,'bIn']), na.rm=TRUE) # p 2397, eq 5
HP <- -rowSums(vlog2v(b/np[,'bOut']), na.rm=TRUE) # p 2397, eq 6
# See comment at bottom of 2397 about about top-level threshold
fracT.q.prime <- mean(np[,'d.prime']>=top.level.threshold)
fracI.q.prime <- mean(0<np[,'d.prime'] & np[,'d.prime']<top.level.threshold)
fracB.q.prime <- mean(0==np[,'d.prime'])
fracT.q <- mean(np[,'d']>=top.level.threshold)
fracI.q <- mean(0<np[,'d'] & np[,'d']<top.level.threshold)
fracB.q <- mean(0==np[,'d'])
# Ratios of resources to consumers - p 2398, eq 11 and 12
NP.q.prime <- 2^(-sum(vlog2v(np[,'nP']/sum(np[,'nP'])), na.rm=TRUE)) /
2^(-sum(vlog2v(np[,'nN']/sum(np[,'nN'])), na.rm=TRUE))
NP.q <- 2^(-sum(vlog2v((np[,'bOut']*np[,'nP'])/sum(np[,'bOut']*np[,'nP'])), na.rm=TRUE)) /
2^(-sum(vlog2v((np[,'bIn']*np[,'nN'])/sum(np[,'bIn']*np[,'nN'])), na.rm=TRUE))
# Link properties
# p 2398, eq 13
LD.q.prime <- (sum(np[,'nP']) + sum(np[,'nN'])) / (2*NumberOfNodes(community))
# p 2398, eq 14
LD.q <- (sum(np[,'bOut']*np[,'nP']/sumb, na.rm=TRUE) +
sum(np[,'bIn']*np[,'nN']/sumb, na.rm=TRUE))/2
# p 2398, col 2
C.q.prime <- LD.q.prime/NumberOfNodes(community)
C.q <- LD.q/NumberOfNodes(community)
# ...the sum of the diversity of outflows weighted by the total outflows,
# of the diversity of inflows weighted by the total inflows. Phi can be
# thought of as the average amount of choice in trophic pathways (Ulanowicz
# and Wolff 1991)
Phi <- sum(HP*np[,'bOut']/sumb) + sum(HN*np[,'bIn']/sumb) # p 2398, eq 16
m <- 2^(Phi/2) # Effective connectance per node p 2398, eq 15
# Page 2398, eq 18
PhiAB <- function(A, B)
{
# A and B should be functions that take a community and return node
# names or indices.
A <- A(community)
B <- B(community)
bOut <- rowSums(b)[A]
bIn <- colSums(b)[B]
return (sum((bOut/sumb) * -vlog2v(b[A,B]/bOut), na.rm=TRUE) +
sum((bIn/sumb) * -vlog2v(t(b[A,B])/bIn), na.rm=TRUE))
}
fracTI.q <- PhiAB(IntermediateNodes, TopLevelNodes) / Phi
fracTB.q <- PhiAB(BasalNodes, TopLevelNodes) / Phi
fracII.q <- PhiAB(IntermediateNodes, IntermediateNodes) / Phi
fracIB.q <- PhiAB(BasalNodes, IntermediateNodes) / Phi
# Page 2399, top left
PhiPrime <- sum(HP/NumberOfNodes(community)) + sum(HN/NumberOfNodes(community))
PhiABPrime <- function(A, B)
{
# A and B should be functions that take a community and return node
# names or indices.
A <- A(community)
B <- B(community)
bOut <- rowSums(b)[A]
bIn <- colSums(b)[B]
s <- NumberOfNodes(community)
return (sum((1/s) * -vlog2v(b[A,B]/bOut), na.rm=TRUE) +
sum((1/s) * -vlog2v(t(b[A,B])/bIn), na.rm=TRUE))
}
fracTI.q.prime <- PhiABPrime(IntermediateNodes, TopLevelNodes) / PhiPrime
fracTB.q.prime <- PhiABPrime(BasalNodes, TopLevelNodes) / PhiPrime
fracII.q.prime <- PhiABPrime(IntermediateNodes, IntermediateNodes) / PhiPrime
fracIB.q.prime <- PhiABPrime(BasalNodes, IntermediateNodes) / PhiPrime
# Weighted chain lengths, p 2399, eq 19 and 20
tlps <- TLPS(community, link.properties=weight)
chains <- TrophicChains(community)
# The biomass flowing through every link in every chain
bc <- apply(as.matrix(chains), 1, function(r)
{
links <- r[which(""!= r)]
sapply(2:length(links), function(index)
{
tlps[tlps$resource==links[index-1] & tlps$consumer==links[index],3]
})
})
cl.q.prime <- 2^(-sapply(bc, function(b.chain) sum(vlog2v(b.chain/sum(b.chain)))))
bc.mean <- sapply(bc, sum)/cl.q.prime
cl.q <- cl.q.prime*bc.mean*nrow(chains) / sum(bc.mean)
# Quantitative generality and vulnerability, p 2400, eq 24-27
nT <- length(TopLevelNodes(community))
nI <- length(IntermediateNodes(community))
nB <- length(BasalNodes(community))
G.q.prime <- sum(np[,'nN'])/(nT+nI)
G.q <- sum(np[,'nN']*np[,'bIn']/sumb, na.rm=TRUE)
V.q.prime <- sum(np[,'nP'])/(nI+nB)
V.q <- sum(np[,'nP']*np[,'bOut']/sumb, na.rm=TRUE)
# Bersier et al table 2
tlps <- TLPS(community, node.properties=c('IsTopLevelNode', 'IsIntermediateNode', 'IsBasalNode'))
fracTI <- with(tlps, sum(resource.IsIntermediateNode & consumer.IsTopLevelNode))/nrow(tlps)
fracTB <- with(tlps, sum(resource.IsBasalNode & consumer.IsTopLevelNode))/nrow(tlps)
fracII <- with(tlps, sum(resource.IsIntermediateNode & consumer.IsIntermediateNode))/nrow(tlps)
fracIB <- with(tlps, sum(resource.IsBasalNode & consumer.IsIntermediateNode))/nrow(tlps)
# This is far faster than ChainLengths
chains.length <- TrophicChainsStats(community)$chain.lengths
Qualitative <- c(FractionTopLevelNodes(community),
FractionIntermediateNodes(community),
FractionBasalNodes(community),
sum(NumberOfConsumers(community)>0) / sum(NumberOfResources(community)>0),
LinkageDensity(community),
DirectedConnectance(community),
fracTI,
fracTB,
fracII,
fracIB,
mean(chains.length),
median(chains.length),
sd(chains.length),
max(chains.length),
Omnivory(community, level=ChainAveragedTrophicLevel),
mean(TrophicGenerality(community)[NumberOfResources(community)>0]),
mean(TrophicVulnerability(community)[NumberOfConsumers(community)>0]),
sd(NormalisedTrophicGenerality(community)),
sd(NormalisedTrophicVulnerability(community)))
Unweighted <- c(fracT.q.prime,
fracI.q.prime,
fracB.q.prime,
NP.q.prime,
LD.q.prime,
C.q.prime,
fracTI.q.prime,
fracTB.q.prime,
fracII.q.prime,
fracIB.q.prime,
mean(cl.q.prime),
median(cl.q.prime),
sd(cl.q.prime),
max(cl.q.prime),
mean(np[,'o.prime']),
G.q.prime,
V.q.prime,
sd(np[,'g.prime']),
sd(np[,'v.prime']))
Weighted <- c(fracT.q,
fracI.q,
fracB.q,
NP.q,
LD.q,
C.q,
fracTI.q,
fracTB.q,
fracII.q,
fracIB.q,
mean(cl.q),
median(cl.q),
sd(cl.q),
max(cl.q),
mean(np[,'o']),
G.q,
V.q,
sd(np[,'g']),
sd(np[,'v']))
res <- cbind(Qualitative, Unweighted, Weighted)
rownames(res) <- c('Fraction top level',
'Fraction intermediate', 'Fraction basal', 'Ratio resources:consumers',
'Link density', 'Connectance', 'Fraction links top:intermediate',
'Fraction links top:basal', 'Fraction links intermediate:intermediate',
'Fraction links intermediate:basal', 'Mean chain length',
'Median chain length', 'SD chain length', 'Max chain length',
'Degree of omnivory', 'Generality', 'Vulnerability',
'SD standardised generality', 'SD standardised vulnerability')
return (res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.