Bipartite network metrics from scratch

Connectance, generality, vulnerability and the specialisation indices H2’ and d’ for a plant-pollinator matrix, all coded in base R, with the size and sampling caveats spelled out.

A plant-pollinator dataset, a host-parasite rearing table or a seed-disperser log all share a shape: a rectangular matrix where rows are one set of species, columns are another, and each cell counts the interactions between a pair. A handful of numbers are then meant to summarise the whole web: how connected it is, how specialised the species are, whether a few players carry most of the links. This post codes those descriptors by hand, so the definitions stay visible, and it flags the two properties that make raw network numbers hard to compare: they depend on the size of the matrix and on how hard the web was sampled.

The reference implementation for all of this is the bipartite package (Dormann et al 2009); here the render path is base R only, which keeps every index a few lines of readable code.

A synthetic plant-pollinator web

The matrix comes from a simple model: each pollinator has an abundance and a niche position, each plant a position on the same axis, and the expected number of visits falls off with the distance between them. Poisson counts finish it off.

set.seed(4260)
A0 <- 12L; P0 <- 10L
poll_ab   <- exp(rnorm(A0, 0, 0.7))
plant_ab  <- exp(rnorm(P0, 0, 0.6))
poll_mu   <- runif(A0, 0, 1)
plant_pos <- seq(0.05, 0.95, length.out = P0)
niche_w   <- 0.18
lam <- outer(poll_ab, plant_ab) * exp(-(outer(poll_mu, plant_pos, "-")^2) / (2 * niche_w^2))
lam <- lam / sum(lam) * 600
W <- matrix(rpois(A0 * P0, lam), A0, P0)
rownames(W) <- paste0("poll", seq_len(A0)); colnames(W) <- paste0("plant", seq_len(P0))
W <- W[rowSums(W) > 0, colSums(W) > 0, drop = FALSE]
A <- nrow(W); P <- ncol(W)
m <- sum(W); rs <- rowSums(W); cs <- colSums(W)
c(pollinators = A, plants = P, interactions = m)

 pollinators       plants interactions 
          12           10          643 

The web has 12 pollinators, 10 plants and 643 recorded interactions.

Connectance and links: the descriptors that scale with size

The most basic number is connectance: the fraction of possible links that are realised. With A rows and P columns there are A * P possible cells, and connectance is the count of non-empty cells divided by that.

n_links   <- sum(W > 0)
connect   <- n_links / (A * P)
links_per <- n_links / (A + P)
c(links = n_links, connectance = round(connect, 3), links_per_species = round(links_per, 2))

            links       connectance links_per_species 
           69.000             0.575             3.140 

Here 69 of the 120 cells are filled, so connectance is 0.575 and there are 3.14 links per species. These are honest counts, but they are the wrong things to compare across studies. Connectance falls as a web grows, because the number of possible cells grows faster than the number a species can plausibly use, and it climbs when a web is sampled harder, because more rare links get recorded. Two webs with identical ecology but different dimensions or effort will report different connectance. The indices below are built to lean on this less.

Quantitative descriptors: effective partners

A pollinator that splits its visits evenly across five plants is more of a generalist than one that makes ninety per cent of its visits to a single plant, even if both touch five plants. The effective number of partners captures that by exponentiating the Shannon diversity of a species’ interactions.

shannon_eff <- function(x) { p <- x[x > 0] / sum(x); exp(-sum(p * log(p))) }
eff_plants <- apply(W, 1, shannon_eff)   # effective plants per pollinator
eff_polls  <- apply(W, 2, shannon_eff)   # effective pollinators per plant
generality   <- sum((rs / m) * eff_plants)
vulnerability<- sum((cs / m) * eff_polls)
link_density <- (generality + vulnerability) / 2
w_connect    <- link_density / (A + P)
round(c(generality, vulnerability, link_density, w_connect), 3)

[1] 4.968 5.017 4.992 0.227

Averaging these effective counts across species, weighted by how many interactions each carries, gives the quantitative descriptors of Bersier et al (2002): generality 4.97 (the effective plants an average pollinator uses), vulnerability 5.02 (the effective pollinators an average plant receives), their mean the weighted linkage density 4.99, and the weighted connectance 0.227. Because they count effective rather than raw partners, a single stray visit barely moves them, unlike the binary link count.

Network specialisation: H2’

The index that travels best across webs is H2’ (Bluthgen et al 2006). It asks how far the observed spread of interactions departs from what the row and column totals alone would produce. Write every cell as a proportion of the grand total and take the two-dimensional Shannon entropy H2. The most even web consistent with the same marginal totals is the one where interactions are independent of identity, and its entropy is just the sum of the marginal entropies; the least even is the most tightly packed web with those totals. H2’ rescales the observed entropy between the two.

pmat  <- W[W > 0] / m
H2    <- -sum(pmat * log(pmat))
r <- rs / m; cc <- cs / m
H2max <- -sum(r[r > 0] * log(r[r > 0])) - sum(cc[cc > 0] * log(cc[cc > 0]))
# least-even coupling with the same marginals: greedily pack the largest feasible cell
pack_entropy <- function(rw, cl, tot) {
  ent <- 0
  repeat {
    if (sum(rw) <= 0) break
    best <- -1; bi <- bj <- NA
    for (i in which(rw > 0)) for (j in which(cl > 0)) {
      v <- min(rw[i], cl[j]); if (v > best) { best <- v; bi <- i; bj <- j }
    }
    p <- best / tot; ent <- ent - p * log(p)
    rw[bi] <- rw[bi] - best; cl[bj] <- cl[bj] - best
  }
  ent
}
H2min   <- pack_entropy(rs, cs, m)
H2prime <- (H2max - H2) / (H2max - H2min)
round(c(H2 = H2, H2max = H2max, H2min = H2min, H2prime = H2prime), 3)

     H2   H2max   H2min H2prime 
  3.741   4.304   2.440   0.302 

The observed entropy 3.74 sits between a maximum of 4.30 and a minimum of 2.44, giving H2’ = 0.30: moderate specialisation. The value runs from 0 for a web no more structured than its totals demand, to 1 for a web where each species interacts only with its own private partners.

To see that H2’ tracks specialisation rather than size, spread the same marginal totals into their independence coupling, which is the most even web those totals allow, and recompute.

W_even <- outer(rs, cs) / m
pe <- W_even[W_even > 0] / sum(W_even)
H2_even <- -sum(pe * log(pe))
H2prime_even <- (H2max - H2_even) / (H2max - H2min)
round(c(H2prime_observed = H2prime, H2prime_evened = H2prime_even), 3)

H2prime_observed   H2prime_evened 
           0.302            0.000 

The evened web returns H2’ = 0.00, exactly zero by construction, while the observed web keeps its 0.30. The two webs share every row and column total; only the internal spread differs, and H2’ reads that difference.

Species specialisation: d’

The single-species version is d’, the standardised departure of one species’ interactions from the availability of its partners. Availability is the column total of each plant as a share of all interactions; a pollinator that tracks availability is a generalist, one that piles onto a rare plant is a specialist. The raw departure is a Kullback-Leibler divergence, rescaled against the most specialised pattern the same species could show.

avail <- cs / m
d_raw <- function(a) { p <- a[a > 0] / sum(a); q <- avail[a > 0]; sum(p * log(p / q)) }
d_max <- -log(min(avail))               # concentrate all interactions on the rarest plant
dprime_poll <- apply(W, 1, function(a) d_raw(a) / d_max)
round(range(dprime_poll), 3); round(mean(dprime_poll), 3)

[1] 0.061 0.387

[1] 0.216

Across pollinators d’ runs from 0.06 to 0.39, averaging 0.22. A value near zero means the species visits plants roughly in proportion to how common those plants are; a value near one means it concentrates on partners far rarer than chance would give.

Figure 1: Left: the interaction matrix, shaded by number of visits. Right: the effective number of plant partners for each pollinator, ranked; the network generality is their interaction-weighted mean.

The honest limit: sampling depth leaks into specialisation

The reason H2’ and d’ are preferred over connectance is that they standardise against the marginal totals, which soaks up part of the abundance and effort signal. Part, not all. A species seen only a few times cannot look like a generalist: with two or three interactions its d’ is pulled upward whether or not it is genuinely choosy. In this web the effect is plain.

conf <- cor(log(rs), dprime_poll)
round(conf, 3)

[1] -0.875

The correlation between a pollinator’s log total interactions and its d’ is -0.87: the pollinators with the fewest recorded visits are the ones scored as most specialised. Some of that is real, because narrow-niche pollinators do interact less in this model, but from the matrix alone the genuine specialists and the under-sampled species cannot be told apart, and both inflate d’. This is the recurring caution in network analysis (Bluthgen 2010; Frund et al 2016), and it sets up the next two posts.

`geom_smooth()` using formula = 'y ~ x'

Figure 2: Left: d’ for each pollinator, ranked and coloured by effective partners. Right: apparent d’ against how many interactions the pollinator was seen making; poorly sampled species trend toward higher specialisation.

Where this leaves you

Connectance and link counts describe a web but cannot be compared across webs of different size or effort. The quantitative descriptors, generality and vulnerability, count effective partners and shrug off stray links. H2’ and d’ standardise against the marginal totals and travel better, but they are still not immune to sampling: a thinly seen species drifts toward apparent specialisation. Two questions follow. Is an observed value larger than a random web with the same totals would give, which is the job of a null model? And how much does the number move as sampling changes? Those are the subjects of the next posts on null models for interaction networks and checking a network analysis.

References

Bersier, L.-F., Banasek-Richter, C. and Cattin, M.-F. (2002). Quantitative descriptors of food-web matrices. Ecology, 83(9), 2394-2407.

Bluthgen, N., Menzel, F. and Bluthgen, N. (2006). Measuring specialization in species interaction networks. BMC Ecology, 6, 9. https://doi.org/10.1186/1472-6785-6-9

Bluthgen, N. (2010). Why network analysis is often disconnected from community ecology: a critique and an ecologist’s guide. Basic and Applied Ecology, 11(3), 185-195. https://doi.org/10.1016/j.baae.2010.01.001

Dormann, C. F., Frund, J., Bluthgen, N. and Gruber, B. (2009). Indices, graphs and null models: analyzing bipartite ecological networks. The Open Ecology Journal, 2, 7-24. https://doi.org/10.2174/1874213000902010007

Frund, J., McCann, K. S. and Williams, N. M. (2016). Sampling bias is a challenge for quantifying specialization and network structure: lessons from a quantitative niche model. Oikos, 125(4), 502-513. https://doi.org/10.1111/oik.02256

Related tutorials

• Nestedness and the NODF metric
• Co-occurrence null models
• Modularity in ecological networks
• Null models for interaction networks