Trophic position/height in food webs (following paths in networks)

Question

As part of developing a demo for a package I'm working on, I need to quantify a classic ecological food web as described below. I have checked out vegan, bipartite and sna but don't see anything that does what I need, though I could be wrong - those are big packages. Thus I'm wondering if this idea is already in package, or if someone has a clever way to compute the results. Seems like it should be a package.

A food web can be described by a matrix of interactions between species A:F, as shown in the code and diagram. In words, one would say that "A eats B which eats E" etc (very hard to see in the matrix, trivial in the diagram).

species <- LETTERS[1:6]

links <- c(0, 1, 0, 0, 0, 0,
    1, 0, 1, 1, 1, 0,
    0, 1, 0, 0, 1, 0,
    0, 1, 0, 0, 1, 1,
    0, 1, 1, 1, 0, 0,
    0, 0, 0, 1, 0, 0)

L <- matrix(links, nrow = 6, byrow = TRUE,
    dimnames = list(species, species))

I want to compute the trophic position and the trophic height of each species. Trophic position is defined as the total number of species the in the food chain below a particular species + 1. In the diagram, the trophic position of A is 6 and for D it is 3. On the other hand, trophic height is the average of the position of a species in each separate chain in which it participates. Species B is connected to 4 different chains (paths); it's height is the average of the positions considered one at at time: (3 + 3 + 3 + 2)/4 = 2.75.

Computationally, one needs to read the matrix L, then work through the different paths implied by the matrix to compute the needed values.

If this is not too obtuse, does anyone know a package that will do this, or see a way to follow the paths and compute on the various lengths/options? It "feels" like there must be some recursive/apply approach that should work, but I don't want to reinvent stuff.

Thanks in Advance

Answer 1

This is classic graph theory stuff. What you have is a directed graph, where animals are nodes and 'A eats B' is an edge. Your matrix is then an adjacency matrix encoding the edges (but you probably want half of the matrix to be 0s, so you don't have A eats B and B eats A. Similarly, zeroes on the diagonal unless you have cannibalistic behaviour). This is a directed acyclic graph as linked in another answer.

Packages for graph theory are listed in the graphical models Task View:

http://ftp.heanet.ie/mirrors/cran.r-project.org/web/views/gR.html

and I think igraph will do what you want. First make a graph object from your matrix (using the modified matrix with only the upper triangle):

> LG=graph.adjacency(L)
> LG
Vertices: 6 
Edges: 7 
Directed: TRUE 
Edges:

[0] 'A' -> 'B'
[1] 'B' -> 'C'
[2] 'B' -> 'D'
[3] 'B' -> 'E'
[4] 'C' -> 'E'
[5] 'D' -> 'E'
[6] 'D' -> 'F'

Then you can use neighborhood.size to get the trophic position:

> neighborhood.size(LG,Inf,mode="out")
[1] 6 5 2 3 1 1

which should be in order A,B,C,D,E,F. Note use of Inf to stop limiting the neighbourhood, and the use of mode="out" to only follow the graph edges in the direction of the eating.

All this talk of food chains is making me hungry. Hope that gets you started.

Barry

Answer 2

Here is one approach to calculating the trophic heights you're after.

Let's say that matrix L encodes a DAG with links FROM the species in each row TO the species in each column. L_(i,j)=1 indicates that spp_i eats spp_j.

Then L*L expresses the number of two-step 'trophic paths' connecting each pair of species, L*L*L contains the number of three step paths, and so on. Between them, the set of 'matrix powers' records all paths, of whatever length, connecting pairs of nodes.

From your description, a species' trophic height is one plus the average path length of all the paths that connect it to one of the 'leaf nodes' at the tips of the graph.

## Here I've edited the matrix L to make it a DAG
species <- LETTERS[1:6]

links <-
  c(0, 1, 0, 0, 0, 0,
    0, 0, 1, 1, 1, 0,
    0, 0, 0, 0, 1, 0,
    0, 0, 0, 0, 1, 1,
    0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0)

L <- matrix(links, nrow = 6, byrow = TRUE,
    dimnames = list(FROM=species, TO=species))

This function (which uses an operator in the expm package) should do what you are looking for. It could do with some more detailed commenting, but I'll add that later if I have the time.

library(expm) ## provides "%^%", a matrix power operator

calcHeight <- function(MAT) {
    ## Find 'leaf nodes' (i.e. species that are only eaten,
    ## and don't eat any others)
    leaves <- which(rowSums(L)==0)

    ## Find the maximum possible chain length (if this is a DAG)
    maxHeight <- nrow(MAT) - length(leaves) - 1
    ## Then use it to determine which matrix powers we'll need to calculate.
    index <- seq_len(maxHeight)

    paths <- lapply(index, FUN=function(steps) MAT %^% steps)
    pathSteps <- lapply(index, FUN=function(steps) (1 + steps) * paths[[steps]])

    ## Trophic height is expressed relative to leaf nodes
    paths <- Reduce("+", paths)[-leaves, leaves]
    pathSteps <- Reduce("+", pathSteps)[-leaves, leaves]

    rowSums(pathSteps)/rowSums(paths)
}

calcHeight(L)
   A    B    C    D 
3.75 2.75 2.00 2.00