Bipartite
This was been updated Friday afternoon, 17 Nov. 2017
# for plotting
library(ggplot2)
# for our web construction and food web properties
library(bipartite)
library(igraph)
library(rARPACK)
# for structural equation models (including path analysis)
library(lavaan)Comparing Consumer vs. Mutualistic Bipartite Webs
In this lab, we’ll take a stab at investiagting the same questions as in Thebault and Fontaine (2010). Are consumptive and mutualistic networks structured differently? If so, how do those structures contribute to the stability of those webs?
To investigate those questions, we will use a different approach than that of Thebault and Fontaine. Thebault and Fontaine (2010) simulated dynamic communities using systems of ODEs that are similar to those that we have seen in earlier in this course. For instance, they used the same sort of mathematical expressions, such as negative density dependence and type II functional responses, that we saw previously in our lab using Rosenzweig and MacArthur’s predator-prey model, and in mutualisms.
We will use the approach of May (1973) and Pimm and Lawton (1977, 1978). In this approach, we do not simulate the dynamics. Instead, we assume that the species involved coexist, and that we know how they interact. Specifically, that we know their instantaneous per capita effect on each other’s population growth rate. That is, we know the Jacobian matrix that describes their interactions. That is what Stuart Pimm and John Lawton assumed in their approach in their groundbreaking research (Stevens, 2009).
In our approach, we will establish interactions at random, in which each pair of species (plant and animal) has the same probability of interacting. That probability is connectance (\(L/L^2\)). This means that the structure (who is connected to whom) will be determined randomly rather than as a result of the dynamics of ODE models as in Thebault and Fontaine (2010). Further, the number of species will be and the connectance will both be draw at random, from specified ranges.
Creating bipartite webs
We’ll start by creating a bipartite web of mutualists: plants and their pollinators.
# The number of species of plants
Sp <- 5
# The number of species of pollinators
Sa <- 12
# the total number of species
S <- Sa+SpWe then set the connectance, and then generate random interactions between species.
# The connectance in this web.
contc <- .5
# Create interactions at random, drawing 0's and 1's with a probability of contc=0.5
interactions <- rbinom( n = Sa * Sp, size=1, prob=contc)
# Put those connections into our web
bn <- matrix(interactions, nr=Sp, nc=Sa)
rownames(bn) <- paste("plant", LETTERS[1:Sp], sep="")
colnames(bn) <- paste("bug", LETTERS[1:Sa], sep="")
#...and look at it
bn ## bugA bugB bugC bugD bugE bugF bugG bugH bugI bugJ bugK bugL
## plantA 0 1 0 1 0 0 1 0 0 0 1 0
## plantB 0 0 0 1 0 0 1 1 1 0 1 1
## plantC 0 0 0 1 0 0 1 0 1 1 0 0
## plantD 1 1 1 1 1 0 1 0 1 0 0 0
## plantE 0 1 0 0 1 0 0 0 1 1 0 1Now that we have randomly constructed web, let’s create a realistic assumption about how much each of the pollinators interact with those plant species. It is generally observed that most interactions are weak (i.e. infrequent) and a few are strong (i.e. frequent). That pattern is an exponential distribution. We will create those interactions and then plot them.
# Exponentially distributed of interaction strengths
# Random numbers from and exponential distribution
x <- rexp(Sa*Sp)
# express these as a probability distribution (divide by the sum)
xp <- x/sum(x)
# assign those exponential probabilites to our interactions.
bne <- bn * xp
# graph everything...
{
layout(matrix(c(1,2,3,3), nc=2, byrow=T) )
hist(bne, main="", xlab="interaction strength")
visweb(bne)
plotweb(bne)
}
That is just an example of a web we will create in our simulation—we create a bipartite web with a random number of species, and random connectance, each within specified bounds.
Investigating our ideas
We have a function we can use to investigate our ideas, which we load here. Make sure that bip_stability.R is in your working directory.
source("bip_stability.R") Simulating many webs
In our webs, we can compare the effects of negative vs. positive interactions directly, because we will create one network with interaction strengths, and then simply reverse the signs of the effects of the animals on the plants to get herbivory vs. mutualism.
We simulate 1000 random bipartite webs using a large number of different numbers of plant and animal species. We select a connectance by drawing randomly from a uniform distribution between the values in argument C.
The function bip_stability() performs eigananalysis on the resulting Jacobian matrix to estimate resilience (the real part of the dominant eigenvalue x -1). It also calculates nestedness and modularity in ways similar to that in Thebault and Fontaine (2010), and it tallies diversity as the total number of species and the connectance used for each of the simulated webs.
b <- bip_stability(Spp.a = 16:60, Spp.p = 8:30, C=c(.05, 0.5), reps=1000)We can get a summary of the output of the 1000 simulated webs here.
summary(b$out)## S C modularity nestedness
## Min. :18.00 Min. :0.06861 Min. :0.0000 Min. : 1.681
## 1st Qu.:43.00 1st Qu.:0.17644 1st Qu.:0.0000 1st Qu.: 9.161
## Median :55.00 Median :0.27598 Median :0.0885 Median :14.819
## Mean :55.38 Mean :0.28297 Mean :0.1930 Mean :14.588
## 3rd Qu.:68.00 3rd Qu.:0.38889 3rd Qu.:0.3644 3rd Qu.:20.296
## Max. :90.00 Max. :0.54545 Max. :0.8289 Max. :30.228
## resilience.m resilience.h
## Min. :0.00919 Min. :0.004902
## 1st Qu.:0.02857 1st Qu.:0.014018
## Median :0.03895 Median :0.018856
## Mean :0.04443 Mean :0.021689
## 3rd Qu.:0.05582 3rd Qu.:0.027152
## Max. :0.17657 Max. :0.077294Are these webs stable? What are the medians and ranges of each of the web properties? How do the resiliences of the two interaction types differ?
We can also graph these relations.
pairs(b$out)
Do any of these relate to figures in Thebault and Fontaine (2010)?
Path Analysis
Now we do the same sort of analysis that Thebault and Fontaine did to get their Fig. 2 C,D.
We start by scaling our data so that all variables have a mean of zero and a standard deviation of 1.
dat<- data.frame( scale(b$out) )We specify each of the relationships, where y ~ x is “y as a function of x”. We then fit the model and get the path coefficients.
We start with the mutualism web.
# Step 1: Specify model
model <- 'resilience.m ~ nestedness + modularity + S + C
nestedness ~ S + C
modularity ~ S + C'
# Step 2: Estimate model of mutualism
mod.m.fit <- sem(model, data=dat)
summary(mod.m.fit)## lavaan 0.6-3 ended normally after 40 iterations
##
## Optimization method NLMINB
## Number of free parameters 11
##
## Number of observations 1000
##
## Estimator ML
## Model Fit Test Statistic 111.512
## Degrees of freedom 1
## P-value (Chi-square) 0.000
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## resilience.m ~
## nestedness 0.297 0.080 3.716 0.000
## modularity 0.178 0.040 4.501 0.000
## S -0.621 0.023 -27.444 0.000
## C 0.420 0.085 4.930 0.000
## nestedness ~
## S 0.032 0.007 4.285 0.000
## C 0.967 0.007 130.177 0.000
## modularity ~
## S -0.314 0.015 -20.930 0.000
## C -0.777 0.015 -51.791 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .resilience.m 0.343 0.015 22.361 0.000
## .nestedness 0.054 0.002 22.361 0.000
## .modularity 0.219 0.010 22.361 0.000Look at the estimates of the regressions. Those are the path coefficients.
Now it is up to you to get a pencil and paper and draw the picture….
Next, do the same thing for herbivory.
# Step 1: Specify model
model <- 'resilience.h ~ nestedness + modularity + S + C
nestedness ~ S + C
modularity ~ S + C'
# Step 2: Estimate model of herbivory
mod.h.fit <- sem(model, data=dat)
hfit <- summary(mod.h.fit) ## lavaan 0.6-3 ended normally after 35 iterations
##
## Optimization method NLMINB
## Number of free parameters 11
##
## Number of observations 1000
##
## Estimator ML
## Model Fit Test Statistic 111.512
## Degrees of freedom 1
## P-value (Chi-square) 0.000
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## resilience.h ~
## nestedness 0.501 0.069 7.222 0.000
## modularity -0.087 0.034 -2.533 0.011
## S -0.780 0.020 -39.618 0.000
## C 0.005 0.074 0.068 0.946
## nestedness ~
## S 0.032 0.007 4.285 0.000
## C 0.967 0.007 130.177 0.000
## modularity ~
## S -0.314 0.015 -20.930 0.000
## C -0.777 0.015 -51.791 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .resilience.h 0.259 0.012 22.361 0.000
## .nestedness 0.054 0.002 22.361 0.000
## .modularity 0.219 0.010 22.361 0.000