library("igraph") # To handle the webs objects and plot them
library("magrittr") # Provides pipes ( %>% )
library("NetIndices") # To compute ecological indices
library("dplyr") # Filter, select, mutate, ... data sets
library("reshape2") # To organise and 'shape' our data sets
library("ggplot2") # Nice plots

options(stringsAsFactors = FALSE) # Ask me why

Ryan F. Hechinger, Kevin D. Lafferty, John P. McLaughlin, Brian L. Fredensborg, Todd C. Huspeni, Julio Lorda, Parwant K. Sandhu, Jenny C. Shaw, Mark E. Torchin, Kathleen L. Whitney, and Armand M. Kuris (2011). Food webs including parasites, biomass, body sizes, and life stages for three California/Baja California estuaries. Ecology 92:791—791. http://dx.doi.org/10.1890/10-1383.1

Web.source <- "http://esapubs.org/archive/ecol/E092/066/CSMweb_Links.txt"

Web.source %>%
  url %>% # Open the connection
  # and read the data
  read.csv(sep = "\t") -> Raw.web.data

Raw.web.data %>% names

##  [1] "PresentAtCSM"              "PresentAtEPB"             
##  [3] "PresentAtBSQ"              "PresentAtAnyEstuary"      
##  [5] "ConsumerNodeID"            "ResourceNodeID"           
##  [7] "ConsumerSpeciesID"         "ResourceSpeciesID"        
##  [9] "ConsumerSpeciesID.StageID" "ResourceSpeciesID.StageID"
## [11] "LinkTypeID"                "LinkType"                 
## [13] "LinkEvidence"              "LinkEvidenceNotes"        
## [15] "LinkFrequency"             "LinkN"                    
## [17] "DietFraction"              "ConsumptionRate"          
## [19] "VectorFrom"                "PreyFrom"

"http://esapubs.org/archive/ecol/E092/066/CSMweb_Links.txt" %>%
  url %>% # Open the connection
  read.csv(sep = "\t") %>% # Read in the raw data
  filter(LinkType %in% "predation") %>% # Filter only the relevant interactions
  filter(LinkEvidence %in% "observed") %>% # and only if directly observed
  select(ResourceNodeID,ConsumerNodeID) %>% # Select the column of From and To
  graph.data.frame(directed=TRUE) -> CSM.graph # And build a igraph web

What do we have?

CSM.graph %>% vcount # How many vertices in the graph?

## [1] 80

CSM.graph %>% ecount # How many edges in the graph?

## [1] 242

Connectance:

Link.density <- function(x) x %>% {ecount(.) / vcount(.)^2}
CSM.graph %>% Link.density

## [1] 0.0378125

Paul Erdős and Alfréd Rényii. (1959)
On Random Graphs. I Publicationes Mathematicae 6: 290–297.

• \(n\) vertices
• each link has i.i.d. probability \(p\)

Hence, each and every graph with \(M\) links appears with probability \[p^M \left( 1 - p \right)^{n^2 - M}\]

CSM.graph %>%
  erdos.renyi.game(n = vcount(.), # Random graph with the same number of nodes...
                   p.or.m = Link.density(.), # ...and link density as the observed
                   type = "gnp", # This is to tell igraph we specify link density
                   directed = T, # We want a directed graph...
                   loops = F # ...with no loops
                  ) -> Random.graph

Shorter version:

CSM.graph %>%
  {erdos.renyi.game(vcount(.), # number of vertices
                    Link.density(.), # link density
                   "gnp", T, F) # type, directed, and no loops
  } -> Random.graph

Notice the curly brackets { } around the game.

For the observed

CSM.graph %>% degree(mode = "total") -> CSM.degree.tot
CSM.graph %>% degree(mode = "in") -> CSM.degree.in
CSM.graph %>% degree(mode = "out") -> CSM.degree.out

And the random

Random.graph %>% degree(mode = "total") -> Random.degree.tot
Random.graph %>% degree(mode = "in") -> Random.degree.in
Random.graph %>% degree(mode = "out") -> Random.degree.out

Prepare the data…

CSM.degrees <- data.frame("graph.name" = "CSM",
                       "degree.tot" = CSM.degree.tot,
                        "degree.in" = CSM.degree.in,
                        "degree.out" = CSM.degree.out)

Random.degrees <- data.frame("graph.name" = "Random",
                          "degree.tot" = Random.degree.tot,
                          "degree.in" = Random.degree.in,
                          "degree.out" = Random.degree.out)

CSM.degrees %>% summary
Random.degrees %>% summary

… massage the data…

CSM.degrees  %<>% melt

## Using graph.name as id variables

Random.degrees %<>% melt

## Using graph.name as id variables

rbind(CSM.degrees,Random.degrees) -> Full.degrees

… define a suitable ggplot function …

Plot.density.degree <- function(data){
  data %>%
    ggplot(aes(x = value, fill=graph.name)) %>%
  + geom_histogram(binwidth = 1, alpha = .5, position = "identity") %>%
  + facet_grid(. ~ variable) %>%
  + theme_bw() %>% return
}

Full.degrees %>% Plot.density.degree

CSM.degree.tot %>% mean

## [1] 6.05

Random.degree.tot %>% mean

## [1] 5.375

Peak.density <- function(array_values) {
  array_values %>% density %$% # compute the density
    x[which.max(y)] %>% return # find the value at which the density is max
}

CSM.degree.tot %>% Peak.density

## [1] 2.776099

Random.degree.tot %>% Peak.density

## [1] 5.015866

Full.degrees %>%
  filter(variable %>% equals("degree.tot")) %>%
  Plot.density.degree %>%
  + geom_vline(xintercept = CSM.degree.tot %>% Peak.density, col = "red") %>%
  + geom_vline(xintercept = Random.degree.tot %>% Peak.density, col = "blue") -> densities.vs.obs

densities.vs.obs

CSM.degree.tot %>% var

## [1] 24.88354

Random.degree.tot %>% var

## [1] 4.870253

\(H_0\): The degree distribution' (mean, peak, variance) is the same under the random (null) model and in the observed food web.

\(H_1\): No, it is not, you silly! How often do you see such a weird distribution!

Idea: sample a random web,
compute its degree distribution,
compute the mean, peak and var,
compare with the observed one,
sample a random web,
compute its degree distribution,
compute the mean, peak and var,
compare with the observed one,
sample a random web,
compute its degree distribution,
compute the mean, peak and var,
compare with the observed one,
sample a random web,
compute its degree distribution,
compute the mean, peak and var,
compare with the observed one,
sample a random web,
compute its degree distribution,
compute the mean, peak and var,
compare with the observed one,
…

Random sampler (we know this already!):

Sampler <- function(Graph) {
  Graph %>% {erdos.renyi.game(vcount(.), # number of vertices
                              Link.density(.), # link density
                              "gnp", T, F) # type, directed, and no loops
             } %>%
  return
}

Compute the degrees:

Degreer <- function(Graph) {
  data.frame("degree.tot" = Graph %>% degree(mode = "total"),
             "degree.in" = Graph %>% degree(mode = "in"),
             "degree.out" = Graph %>% degree(mode = "out")) %>%
  return
}

Compute the desired statistics:

# we define an ancillary function to compute the statistics we want
Stats <- function(Degrees,mode){
  data.frame("Mode" = mode, # We record whether in, out or total
             "Mean" = mode %>% Degrees[,.] %>% mean, # We want to compare the mean
             "Peak" = mode %>% Degrees[,.] %>% Peak.density, # The density peak
             "Var"  = mode %>% Degrees[,.] %>% var) %>% # And the variance
  return
}

For in, out and total degree:

Statistiker <- function(Degree.data) {
  rbind(
    Degree.data %>% Stats("degree.in"),
    Degree.data %>% Stats("degree.out"),
    Degree.data %>% Stats("degree.tot")
  ) %>% # We bind all together
  return
}

Now we have a pipeline from the graph to the statistics:

CSM.graph %>% # Get the graph
  Degreer %>% # Compute the degrees
  Statistiker %>% # Compute the statistics
  melt %>% # Melt the data (because plots)
  # mutate adds a news variable
  mutate(Origin = "observed") -> CSM.stats

Prepare the data…

R <- 42 # How many randomizations?
# Each sample of Statistics() is made by 3 rows
L <- R*3
# We initialise an empty data frame to store the simulations data
data.frame("Mode" = L %>% character,
           "Mean" = L %>% numeric,
           "Peak" = L %>% numeric,
           "Var"  = L %>% numeric) -> Rands.stats

… fill the data.

# Each sample of Statistics() is made by 3 rows
run.range <- seq(from = 1, to = L, by = 3)
for(run in run.range) {
  CSM.graph %>%
    Sampler %>%
    Degreer %>%
    Statistiker -> Rands.stats[run:(run+2),]
}

Transform the data to the same format of CSM.data

Rands.stats %<>%
  melt %>% # Melt the data (because plots)
  # mutate adds a news variable
  mutate(Origin = "random")

Put everything together

Full.stats <- rbind(Rands.stats,CSM.stats)

… define a suitable ggplot function …

Plot.comparison <- function(Data.to.plot){
  Data.to.plot %>%
    ggplot(aes(x = value))%>% # the value on x and density on y
  + geom_histogram( # plot an histogram with the 'random' data
    data = Data.to.plot %>% filter(Origin %>% equals("random")), 
                   binwidth = 1) %>% # and set a proper bin width
  + geom_vline( # plot a vertical line with the 'observed' data
    data = Data.to.plot %>% filter(Origin %>% equals("observed")),
               aes(xintercept = value), # The observed statistics value
               col = "red", size = 2) %>% # And make it thick and red
  + facet_wrap(Mode ~ variable, scale="free") %>% # Do that for each Mode and statistics
  return
}

Full.stats %>% Plot.comparison

• \(x\) axis Omnivory Index (from NetIndices)
• \(y\) axis Trophic Level (from NetIndices)
• vertex size degree (from igraph)
• vertex color centrality (betweenness) rank (from igraph)

Plot.foodweb <- function(Graph, dodge=0.08, heat=10, size.scale=3) {
   if("viridis" %in% rownames(installed.packages()) == FALSE) {install.packages("viridis")}
   library(viridis)
   Layer <- function(Graph){
       Graph %>% vcount -> V # How many vertices in the graph?
       Graph %>% get.adjacency(sparse = F) %>% TrophInd -> Trophs
       # We define a matrix with two columns:
       c(Trophs$OI + V %>% runif(max=dodge), # one with Omnivory Index (and a bit of noise)
         Trophs$TL - 1 # and one with the Trophic Levels
        ) %>% matrix(ncol=2) %>%
       return
     }
  Graph %>% betweenness(di=F,nor=T) %>% rank(ties="min") %>% max %>% viridis -> col_palette
  Graph %>% betweenness(di=F,nor=T) %>% rank(ties="min") %>% col_palette[.] -> V(Graph)$color
  Graph %>% get.adjacency(sparse = F) %>% TrophInd %$% OI -> Omni
  Graph %>% plot.igraph(layout= Layer(.),
                        vertex.label=NA,
                        vertex.size= degree(.) %>% sqrt * size.scale,
                        edge.arrow.size=.5,
                        edge.width=.5)
}

par(mar=c(0,0,0,0))
CSM.graph %>% Plot.foodweb

Everything into a big function!

Tester <- function(Observed.Graph,Runs) {
  L <- Runs *3
  Observed.Graph %>% Degreer %>% Statistiker %>%
    melt %>% mutate(Origin = "observed") -> Observed.stats
  data.frame("Mode" = L %>% character,
             "Mean" = L %>% numeric,
             "Peak" = L %>% numeric,
             "Var"  = L %>% numeric) -> Rands.stats
  run.range <- seq(from = 1, to = L, by = 3)
  for(run in run.range) {
    Observed.Graph %>% Sampler %>%
    Degreer %>% Statistiker -> Rands.stats[run:(run+2),]
  }
  Rands.stats %>% melt %>%
    mutate(Origin = "random") %>% 
    rbind(Observed.stats) %>% Plot.comparison
}

Everything into a big function!

Graph.grabber <- function(web.source,linktype,linkevidence){
  web.source %>%
    url %>% # Open the connection
    read.csv(sep = "\t") %>% # Read in the raw data
    filter(LinkType %in% linktype) %>% # Filter on interactions type
    filter(LinkEvidence %in% linkevidence) %>% # and evidence type
    select(ResourceNodeID,ConsumerNodeID) %>% # Select the column of From and To
    graph.data.frame(directed=TRUE) %>%
return
}

"http://esapubs.org/Archive/ecol/E092/066/EPBweb_Links.txt" %>%
  Graph.grabber(c("predation",
                  "Social Predation",
                  "predation on free-living non-feeding stage"),
                c("observed",
                  "inferred")) %>%
  Plot.foodweb

"http://esapubs.org/Archive/ecol/E092/066/EPBweb_Links.txt" %>%
  Graph.grabber("predation","observed") %>%
  Tester(100)