Spatial Guilds in the Serengeti Food Web Revealed by a Bayesian Group Model

Edward B. Baskerville1,∗, Andy P. Dobson2, Trevor Bedford1,3, Stefano
Allesina4, Mercedes Pascual1,3

1 Department of Ecology and Evolutionary Biology, University of Michigan, Ann
Arbor, Michigan, United States of America
2 Department of Ecology and Evolutionary Biology, Princeton University,
Princeton, New Jersey, United States of America
3 Howard Hughes Medical Institute, University of Michigan, Ann Arbor,
Michigan, United States of America
4 Department of Ecology and Evolution, Computation Institute, The University
of Chicago, Chicago, Illinois, United States of America
$\ast$ E-mail<EMAIL_ADDRESS>
## Abstract

Food webs, networks of feeding relationships among organisms, provide
fundamental insights into mechanisms that determine ecosystem stability and
persistence. Despite long-standing interest in the compartmental structure of
food webs, past network analyses of food webs have been constrained by a
standard definition of compartments, or modules, that requires many links
within compartments and few links between them. Empirical analyses have been
further limited by low-resolution data for primary producers. In this paper,
we present a Bayesian computational method for identifying group structure in
food webs using a flexible definition of a group that can describe both
functional roles and standard compartments. The Serengeti ecosystem provides
an opportunity to examine structure in a newly compiled food web that includes
species-level resolution among plants, allowing us to address whether groups
in the food web correspond to tightly-connected compartments or functional
groups, and whether network structure reflects spatial or trophic
organization, or a combination of the two. We have compiled the major
mammalian and plant components of the Serengeti food web from published
literature, and we infer its group structure using our method. We find that
network structure corresponds to spatially distinct plant groups coupled at
higher trophic levels by groups of herbivores, which are in turn coupled by
carnivore groups. Thus the group structure of the Serengeti web represents a
mixture of trophic guild structure and spatial patterns, in contrast to the
standard compartments typically identified in ecological networks. From data
consisting only of nodes and links, the group structure that emerges supports
recent ideas on spatial coupling and energy channels in ecosystems that have
been proposed as important for persistence. Our Bayesian approach provides a
powerful, flexible framework for the study of network structure; we believe it
will prove instrumental in a variety of biological contexts.

## Introduction

Food webs, networks of feeding relationships in ecosystems, connect the biotic
interactions among organisms with energy flows, thus linking together
population dynamics, ecosystem function, and network topology. Ecologists have
been using this powerful conceptual tool for more than a century [1, 2, 3].
One element of structure of particular relevance to large food webs is the
subdivision of species into compartments or groups, a feature that has been
proposed to contribute to food web stability by constraining the propagation
of disturbances through a network [4]. Although a large literature has
considered the presence and dynamic significance of compartments in food webs,
on the whole, evidence for the importance of compartments has been
inconclusive [5, 6, 7, 8]. In this literature, compartments are alternately
referred to as modules, clusters, or “communities” [9], and are defined by
high link density within groups and low link density between them.

Recent work with a probabilistic model considers a more flexible notion of
groups, allowing link density to be high or low within any group or between
any pair of groups [8]. Groups can thus represent compartments in the previous
sense, but can also represent trophic guilds or roles [10, 11], sets of
species that feed on, and are fed on, by similar sets of species. By fitting
models of this type to data, the dominant topological pattern in the network
can be found, which may include (spatial) compartments, or trophic guilds, or
some combination of the two. The initial application of this model to
empirical food webs from different ecosystems has revealed a predominance of
trophic guilds rather than compartments [8].

Two major challenges limit the application of this model in resolving the
group structure of large food webs and interpreting its biological basis.
First, most food webs have poor resolution of primary producers; plants in
terrestrial systems and phytoplankton in aquatic ones are typically
represented by a few nodes that are highly aggregated taxonomically. This is
unfortunate, because their position at the base of a food web means that
primary producers are essential for understanding how the web is spatially
organized, and how this spatial organization percolates up through the trophic
relationships of other species. Fully resolved primary producers are also
essential in examining how the structure of other trophic levels cuts across
their spatial distribution, and, in so doing, couples different habitats.

Second, some technical problems have hindered the use of probabilistic models
in analyzing group structure. Early food web models served as null models for
food web structure and were tested by generating model webs and comparing
summary statistics against data from real webs [12, 13]. More recently, a more
rigorous approach for measuring the goodness of fit of a model has been
provided by maximum likelihood and model selection [14, 8]; two problems still
remain within this framework. One is technical: standard model-selection
criteria are not applicable to “discrete parameters” such as group membership.
The second problem is more fundamental: there are many almost equally good
arrangements, and it is desirable to extract information not just from a
single best arrangement, but also from the rest of the ensemble.

The Bayesian approach is gaining popularity in ecological modeling due to the
philosophical and conceptual appeal of explicitly considering uncertainty in
parameter estimation as well as its methodological flexibility [15]. This
approach is especially well-suited for handling uncertainty in complex food
web models, and allows us to overcome the limitations of the previous
implementation of the group model. In network inference, there are only a few
examples of complete Bayesian models [16, 17] and a few examples of MCMC for
maximum-likelihood inference [18, 19], but Bayesian inference in phylogenetics
has been long established [20, 21], and provides a clear methodological
analogue.

In this paper, we address the group structure of a newly assembled food web
for the mammals and plants of the Serengeti grassland ecosystem of Tanzania;
we use a new and novel approach to the identification of groups based on
Bayesian inference. We specifically ask whether the structure that emerges
reflects the underlying spatial dimension, as delineated by the different
plant communities that characterize different sub-habitats within the
ecosystem, or whether it is determined by trophic dimensions in the form of
species guilds that share functional roles.

The Serengeti food web is emerging as the most highly resolved terrestrial web
to date [22]. The Serengeti has been studied as an integrated ecosystem for
almost five decades [23, 24, 25], and because of widespread popular
familiarity with the consumer-resource dynamics of lions, hyenas, wildebeest,
zebra and grasses, it provides a strong intuitive test for probabilistic food
web models. Most importantly, at the primary producer level, the Serengeti
food web includes a number of distinct grass and woodland plant communities on
different soils and across a rainfall gradient [26]. The sequential and well-
documented changes in the underlying plant diversity allows us to examine the
extent to which grassland communities define network topology at higher
trophic levels. The high resolution of the species that comprise the mammal
and plant communities of the ecosystem allow us to address the role of space
in the group structure of the food web.

From a complex network of only nodes and links that represent species and
their interactions, the groups that emerge from an otherwise blind
classification of species make remarkable biological sense. We find that the
group structure identified in the Serengeti food web represents a mixture of
trophic guild structure and spatial patterns, an arrangement that differs
significantly from the compartmental structure typically identified in
ecological networks.

## Results/Discussion

### Bayesian Inference and Model Selection for Food Webs

#### Probabilistic Models for Food Webs

In this paper, we use probabilistic modeling as a tool for formalizing
hypotheses about food web structure. We treat a food web, an observed network
of who eats whom in an ecosystem, as data. We start with the basic question:
assuming a probabilistic model of food web structure, what is the probability
of observing this particular real-world food web? This probability is referred
to as the likelihood of observing the data, given model parameters. In a
maximum-likelihood framework, the mechanical part of the inference process is
to find the set of model parameters that makes the likelihood as great as
possible, with the interpretation that this represents the best point estimate
of the underlying process.

We begin with the group model of Allesina and Pascual [8], which was
originally treated in a maximum-likelihood framework. Conceptually, this model
encodes the simple hypothesis that species can be divided into groups, and
species in the same group have statistically similar behavior: they tend to
consume species in certain groups and tend to be consumed by species in
certain groups. Specifically, the probability that a species belonging to
group $i$ is eaten by a species belonging to group $j$ is given by $p_{ij}$,
and conversely, the probability of a link being absent is $(1-p_{ij})$. If
there are $K$ groups, then a matrix $\mathbf{P}$ of $K^{2}$ link probabilities
is required to completely describe the relationships among all groups. The
likelihood for the whole network is the product over all pairs of species of
the probability of a link being present (if present) or absent (if absent). In
the statistical literature, this model structure is known as a stochastic
block model [27]. The assignment of species to groups is also an unobserved
parameter in this model, which adds a layer of difficulty to parameter
estimation. For example, in a network of $100$ species, there are
approximately $5\times 10^{115}$ different ways to partition the network into
groups (see Methods). That is, if you had a computer that could process
$10^{80}$ partitions (as many partitions as there are atoms in the universe)
every femtosecond ($10^{-15}$ s), it would take $1.5\times 10^{13}$ years to
process them all. (By comparison, the universe is only $1.4\times 10^{10}$
years old.)

The group model allows for a more flexible definition of groups than standard
approaches to network clustering, which find groups that have large numbers of
internal connections and relatively few connections between groups [9].
Because each $p_{ij}$ parameter may take any value between $0$ and $1$, good
model fits may result from other relationships, such as high link density
between groups and low link density within groups, and may accommodate
different relationships in different parts of the network. In general, the
best-fitting partitions will try to maximize or minimize the number of links
within specific groups and between specific pairs of groups.

#### Bayesian Inference and Priors for the Group Model

In a Bayesian framework, rules of probability are taken to govern both the
data and model parameters. Rather than finding the set of parameter values
that maximize the likelihood, the goal becomes to estimate a probability
distribution over parameters based on observed data. In this way, we can
directly quantify the uncertainty in our parameters in terms of probabilities.
This permits questions such as: what is the probability that a parameter lies
in a particular range? The name “Bayesian” comes from Bayes’ rule, which tells
us how to use conditional probability statements to infer a posterior
distribution, in this case, the probability distribution over parameter values
conditional on having observed the data, $\mathrm{Pr}(\theta|D)$. If we are
dealing entirely with discrete probability distributions, Bayes’ rule takes
its most intuitive form:

$\displaystyle\mathrm{Pr}(\theta|D)$
$\displaystyle=\frac{\mathrm{Pr}(\theta)\mathrm{Pr}(D|\theta)}{\mathrm{Pr}(D)}.$
(1)

The numerator of the right-hand-side is the probability of producing the data
from the given parameters: the prior probability of those parameters,
$\mathrm{Pr}(\theta)$, times the probability of producing the data given those
parameters, $\mathrm{Pr}(D|\theta)$, the likelihood. The denominator is the
marginal probability of observing the data unconditional on the particular
parameter values at play, which is simply the sum of the probabilities of all
the different ways of producing the data using all possible parameter values,
$\mathrm{Pr}(D)=\sum_{\theta}\mathrm{Pr}(\theta)\mathrm{Pr}(D|\theta)$. In
other words, in order to calculate the posterior probability of parameters
$\theta$, we add up all the different ways of producing the data weighted by
their probability, and then calculate what fraction of that probability came
from parameters $\theta$. From here, we will write these quantities in more
general notation, suitable for a mix of discrete and continuous probability
distributions:

$\displaystyle f(\theta|D)$
$\displaystyle=\frac{f(\theta)f(D|\theta)}{\int_{\theta}f(\theta)f(D|\theta)\,d\theta}$
(2)

where the integral sign represents a multiple integral over discrete and
continuous parameters.

In the Bayesian framework, the model includes not only the formulation of the
likelihood but also a prior distribution over parameters. With the group
model, this means defining a prior distribution over both link probabilities
and arrangements into groups (“partitions”). In general, priors may
incorporate informed knowledge about the system, but in this case we simply
use them to encode different variants of the same basic model. We use two
distributions for partitions and two distributions for link probabilities,
which are combined to form four different model variants.

The two alternative distributions for elements $p_{ij}$ of the link
probability matrix $\mathbf{P}$ are (1) a uniform distribution between 0 and
1, and (2) a beta distribution with shape parameters $\alpha$ and $\beta$,
which are in turn governed by exponential distributions with mean $1$. With
$\alpha$ and $\beta$ fixed at their means, alternative (2) reduces to a
uniform distribution; at other values, the distribution may take a uniform,
convex, concave, or skewed shape. Alternative (2) is thus structured
hierarchically, with exponential hyperpriors for $\alpha$ and $\beta$
governing the beta prior for elements of $\mathbf{P}$.

For partitions, we consider (1) a uniform distribution and (2) a distribution
generated by the Dirichlet process, sometimes referred to as the “Chinese
restaurant process” [28]. Alternative (2) is controlled by an aggregation
parameter $\chi$ that is in turn drawn from an exponential distribution with
mean $1$. The uniform distribution assigns equal prior probability to each
possible partition, irrespective of the number of groups. Because there are
far more ways to partition the network at an intermediate, but relatively
high, number of groups, the uniform prior implicitly biases the model toward
that number. For example, in the Serengeti food web, there are $170$ nodes,
yielding an a priori expectation of $43.8$ groups. In contrast, the
hierarchically structured Dirichlet process prior provides flexibility via the
aggregation parameter $\chi$. When $\chi$ is large, partitions tend to have
many small groups; when $\chi$ is small, partitions tend to have fewer groups,
with a skewed group size distribution.

We also consider two simple models without groups as null comparisons: (1) a
directed random graph model (i.e., one group) with a uniform prior on a single
link probability parameter $p$, and (2) a fully parameterized model, with each
species in its own group (and a $170\times 170$ link probability parameter
matrix $\mathbf{P}$, also with a uniform link probability prior.

For a fuller discussion of models and prior distributions, in particular the
properties of the distribution generated by the Dirichlet process, see
Methods.

#### Bayesian Model Selection via Marginal Likelihood

The Bayesian framework provides a natural way to make probabilistic inferences
based on a particular model. However, we also want to be able to choose
between different models by quantifying their relative goodness of fit. One
approach to Bayesian model selection can be framed directly in terms of Bayes’
rule, mirroring the process for estimating the posterior distribution over
parameters for a single model.

Consider two models, $M_{1}$ and $M_{2}$, to which we assign prior weight
$\mathrm{Pr}(M_{1})$ and $\mathrm{Pr}(M_{2})$. After the data has been
observed, we can calculate the posterior probability of the models using
Bayes’ rule:

$\displaystyle\mathrm{Pr}(M_{1}|D)$
$\displaystyle=\frac{\mathrm{Pr}(M_{1})\mathrm{Pr}(D|M_{1})}{\mathrm{Pr}(D)},$
(3) $\displaystyle\mathrm{Pr}(M_{2}|D)$
$\displaystyle=\frac{\mathrm{Pr}(M_{2})\mathrm{Pr}(D|M_{2})}{\mathrm{Pr}(D)},$
(4)

where the denominator is equal to the probability of observing the data
unconditional of the particular model at play,
$\mathrm{Pr}(D)=\mathrm{Pr}(M_{1})P(D|M_{1})+\mathrm{Pr}(M_{2})P(D|M_{2})$.
The probabilities
$\mathrm{Pr}(D|M_{1})=\int_{\theta_{1}}f(\theta_{1})f(D|\theta_{1})\,d\theta_{1}$
and
$\mathrm{Pr}(D|M_{2})=\int_{\theta_{2}}f(\theta_{2})f(D|\theta_{2})\,d\theta_{2}$
are the marginal likelihoods of the two models, corresponding to the
denominator in Equation 2. If we give the two models equal prior weight, then
the relative posterior weight of the two models is simply given by the
marginal likelihoods. This reasoning extends naturally to any number of
models.

The ratio of the marginal likelihoods is often called the Bayes factor [29,
30, 31], and is equal to the posterior odds ratio of the two models, assuming
equal prior weight:

$\displaystyle B_{12}$
$\displaystyle=\frac{\mathrm{Pr}(D|M_{1})}{\mathrm{Pr}(D|M_{2})}$ (5)

The Bayes factor provides a convenient way to compare models: if $B_{12}=10$,
then we consider support for model $M_{1}$ to be ten times stronger than model
$M_{2}$. In AIC-based model selection, the Bayes factor is analogous to a
ratio of Akaike weights [32].

#### Consensus Partitions

The output of an MCMC simulation includes a long sequence of network
partitions representing draws from the posterior distribution over partitions.
As these partitions are potentially all distinct from each other, but
represent similar tendencies of species to be grouped together, it is useful
to try to summarize the information contained in all the samples in a more
compact form. One approach is to construct a pairwise group-membership matrix
for species in the food web, with entries equal to the posterior probability
that two species are in the same group. A visual representation of this matrix
can illuminate the group structure, and a consensus partition can then be
constructed from this matrix using a simple clustering algorithm. (For more
details, see Methods.)

### Bayesian analysis of the Serengeti food web

#### The Serengeti Data Set

We compiled the food web from published accounts of feeding links in the
literature [33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 26, 43, 44, 45].

The compiled Serengeti food web (Tables S1 and S2 and Figure 1) consists of
$L=667$ feeding links among $S=170$ species ($136$ plants, $25$ herbivores,
and $9$ carnivores). $550$ of the links are herbivorous, and $117$ are
predatory. The fraction of all possible links or connectance ($C=L/S^{2}$),
ignoring all biological constraints, is equal to $0.0231$.

#### Performance of Model Variants

We find unequivocal support for the use of group-based models in describing
the Serengeti food web. Models with group structure have vastly greater
marginal likelihoods than simple null models that ignore group structure
(Table Tables). Furthermore, the use of both the Dirichlet process prior for
network partitions and the beta prior for link probability parameters vastly
improved the fit of the basic model. The best model variant as measured by
marginal likelihood included both the beta prior on link probabilities and the
Dirichlet process prior on partitions. The next best variant included a
uniform partition prior and beta link probability prior, followed by the
variant with Dirichlet process prior and uniform link probability prior. The
strongest variant surpassed its closest competitor by $106$ log-orders in
likelihood and surpassed the model with two uniform priors by $470$ log-orders
in likelihood, providing unequivocal support for including both flexible
priors in the model specification. Accordingly, in the remaining analysis we
consider only the best model variant.

#### Identification of Model Parameters

The posterior mean number of groups $K$ is $14.9$ (95% credible interval
$12,18$), and the mean value of the Dirichlet process parameter $\chi$ is
$3.2$ ($1.6,5.3$) (Figure 2). The prior expectation of $\chi$ was 1.0 and the
prior expectation of $K$ was $5.4$. The finding of posterior values
substantially greater than prior values strongly supports the presence of
detailed group structure in the Serengeti food web.

Mean values for beta distribution parameters are $\alpha=0.046$
($0.029,0.068$) and $\beta=0.89$ ($0.50,1.43$) (Figure S2). The corresponding
beta prior has support concentrated near 0, since most species do not feed on
most other species (Figure S3).

#### Groups Identified in the Serengeti Food Web

The structure of the Serengeti food web is best represented by groups
containing trophically similar species, subdivided by specialization on
different components of prey species, with plant species corresponding to
spatially distinct habitats. The best consensus partition, with $16$ groups,
is shown in Table 2. There are 6 groups of plant species (groups 11–16), 7
groups of herbivore species (groups 4–10) and 3 groups of carnivores species
(groups 1–3). On average, plant groups contain more species than herbivore and
carnivore groups ($22.7$, $3.6$ and $3.0$, respectively). As evident in the
pairwise group-membership matrix (Figure 3), the carnivore and herbivore
groups are well-defined, including several individual species or pairs of
species with distinct diets. Plant groups demonstrate mild overlap, indicating
a partially hierarchical relationship between smaller groups and larger
groups.

Figures 4, 6, and 5 show three alternate views of the food web, organized by
the $16$-group consensus partition. Groups $11$–$13$, $14$–$15$, and $16$
consist of plants located in grassland, woodland, and kopje habitat,
respectively. Group $11$ includes short grasses that are grazed by a wide
variety of species, distinguished from group $12$, which includes short
grasses fed on only by the large grazers ($4$). The 1-species groups (2,
golden jackal Canis aureus; 6, African buffalo Syncerus caffer; 8, agama
lizard Agama planiceps; 10, elephant Loxodenta africana; and 13, Microchloa
kunthii) represent individual species with distinct diets or consumers.

Plant groups are coupled by groups of herbivores, which are in turn coupled by
groups of carnivores. Large migratory grazers ($4$, wildebeest, zebra, and
gazelles) feed primarily on plant groups consisting of grasses ($11$ and
$12$), and are the exclusive consumers of a group of grasses adapted to
seasonal rainfall on sandy or volcanic soils ($12$) that dominate the short-
grass plains in the southern part of the ecosystem. Herbivores feeding in the
longer grasslands and woodlands and in riparian habitats (group $5$) couple
the first grass group ($11$) with woodland plants (group $14$), which are also
consumed to a lesser extent by the large grazers. The hyraxes (group $7$) and
group $9$ (giraffe, olive baboon, and dik-dik) couple kopje habitat (group
$16$) with both woodland and grassland plants. At the highest trophic level,
the large carnivores ($1$) integrate across all the herbivore groups; smaller
carnivores ($2$, $3$) show more specialized diets, reflecting the more
distinct habitats in which they are usually found.

### Discussion

#### Spatial Guilds in the Serengeti Food Web

In order to analyze the group structure of the Serengeti food web, we used a
flexible Bayesian model of network structure that includes no biological
information aside from a set of nodes representing species and links
representing their interactions. The groups that emerge from an otherwise
blind classification of species make remarkable biological sense. Species are
divided into trophic guilds that reveal a clear relationship between the
spatial organization of plant, herbivore, and carnivore groups and the
structure of the network. At the coarsest scale, the groups in the Serengeti
food web correspond to carnivores, herbivores, and plants. The further
subdivisions that emerge within carnivores, herbivores, and plants reveal a
spatial dimension to feeding structure that is only evident because of high
species resolution at the plant level. Ultimately, the group structure we have
derived mirrors the flow of energy up the food web from different spatial
locations, with herbivores integrating spatially separated groups of plants,
and carnivores integrating spatially widespread herbivores. Although the
addition of birds, reptiles, invertebrates, and pathogens will certainly add a
significant number of new groups, we do not expect them to significantly
modify the derived structure for the mammal and plant community. Nor will they
modify the larger tendency for groups to be assembled in ways that reflect the
underlying spatial and trophic structure of the species in the web.

Recently, interesting theoretical and empirical work has highlighted the
relationship between observed patterns of food-web structure and energy flow
that seemingly mirrors the trophic guild structure in the Serengeti. Rooney
and colleagues [46] give evidence that real ecosystems may be dominated by
nested sets of fast and slow “energy channels,” each of which represents a
food chain of trophic guilds. They suggest that this pattern may have a strong
stabilizing effect, based on theoretical work by McCann on spatially coupled
food webs [47]. The group structure for the Serengeti web that emerges from
our analysis supports a pattern of spatial coupling at multiple trophic
levels: the grasslands have very high turnover rates compared to those of the
kopjes and woodlands. This suggests a similar pattern of fast and slow energy
channels to those described by Rooney and colleagues, with fast energy flow up
through the highly seasonal but very productive grasses of the short-grass
plains. These are almost completely consumed by wildebeest and zebra during
their peak calving season, which are then in turn consumed by large predators
(lions and hyenas). In contrast, the resident herbivore species living on
kopjes or in the woodlands reproduce at slower rates and are consumed less
frequently by large carnivores, except during the time when they are unable to
feed on migratory wildebeest and zebra. The group model and the inference
approach presented here allow the examination of the dynamical consequences of
this type of structure to be fully rooted in an empirical pattern,
complementing more theoretical considerations of its central importance for
preserving biodiversity [47].

These patterns emerge directly from the topology of the food web without being
explicitly labeled as different habitats upfront as was done in previous
empirical work [46], showing that topological analysis can reveal structures
that may be very significant for food-web dynamics. They are subtly different,
however, from the proposed pure fast and slow chains, in that they incorporate
the migration of the keystone species in the ecosystem, so the fastest energy
chain is seasonally ephemeral and may only operate for three to four months in
any year. We suspect that even within the sub-habitats of kopjes and woodlands
there are similarly nested faster and slower chains that involve species for
which we are still collating data (e.g., birds, small mammals, and insects).

#### Bayesian Analysis of Food-Web Structure

In this paper, we used a probabilistic model to analyze the structure of a
single food web, an approach we have seen in only one other study based on a
probabilistic version of the niche model [19]. This approach has proved
fruitful in Bayesian phylogenetics, where the combinatorial challenges are
similar. Moreover, we view the group model as only a starting point for richer
modeling efforts to help identify relevant processes that influence the
structure of ecological communities.

In fact, the Bayesian approach described here provides a powerful general
framework for encoding hypotheses about the structure of food webs and
comparing models against each other, and we see it as a natural next step in
the current trend of representing food-web models in a common way. Simple
abstract models such as the niche model and the group model used here act as
proxies for the high-dimensional trait space that determines feeding
relationships in an ecosystem. The identification of actual traits that
correspond to groups (or niche dimensions) is another valuable direction, so
far followed primarily by finding correlations between compartments/groups
[48] or niche values [13] and traits such as body size or phylogenetic
relatedness. Another approach is to directly incorporate these traits into the
probabilistic models, either as model predictors or as informed priors. Both
kinds of analyses are valuable, but the second kind becomes more approachable
in a general Bayesian modeling framework.

The use of flexible hierarchical priors for model parameters is one
straightforward innovation possible in the Bayesian framework. The number of
groups identified by the model increases dramatically with the use of a
flexible beta prior distribution for link probability parameters. In that
model variant, we effectively introduce two degrees of freedom to the model
(the beta distribution parameters) but dramatically reduce the effective
degrees of freedom of the link probability parameters. Note that we properly
penalize parameters by using the marginal likelihood for model selection, so
that the model selection rep- resents a balance between goodness of fit and
model complexity. Moreover, this structure makes intuitive sense: since most
link probability parameters are simply zero, they should not be penalized. An
alternate approach is to remove and add parameters to the model, but this
hierarchical technique is much easier to implement in practice.

Advanced Markov-chain Monte Carlo methods make it possible to accurately
estimate marginal likelihoods for probabilistic network models. Unlike
information criteria such as AIC, BIC, or DIC, an accurate estimate of the
marginal likelihood provides a direct measurement of goodness of fit that
takes into account the degrees of freedom in a model without making any
asymptotic assumptions about parameter distributions [49], and can handle
discrete parameters such as partitioning into groups that are not properly
handled by AIC and BIC.

The Bayesian approach also serves as a means to avoid fundamental issues
inherent in network models with a large parameter space. In a recent study,
Good and colleagues [50] examined the properties of module-finding in networks
using the pervasive modularity-maximization approach [51], finding that even
in relatively small networks a large number of good solutions exist. A
maximization algorithm is thus guaranteed to find a single local maximum of
many—possibly even the best one, but certainly not one that captures the full
range of good solutions. This problem arises whether the quantity to be
maximized is a heuristic such as modularity or a likelihood value. The group
model and other parameter-rich models presumably suffer from similar
degeneracy problems. In the present case, we find that every partition sampled
from the posterior distribution for the best-fitting group model variant is
unique. Although MCMC sampling cannot reproduce the full posterior
distribution, it is an important step in the right direction. Philosophical
arguments aside, one of the main reasons for maximizing likelihood or
modularity is simply that a single solution is far more tractable than a
distribution. The consensus partitioning heuristic used here is one attempt at
recovering a simpler object of study (see Methods); more sophisticated
approaches will be welcomed.

The group model, based on the simple notion that groups of species may have
similar feeding relationships to other groups, reveals that trophic guilds are
the topologically dominant type of group in the Serengeti food web. The model
also reveals an interesting relationship between spatial structure and network
structure that corroborates recent ideas on spatial coupling in food webs. A
theoretical study with a dynamical model suggests that this type of structure
may contribute to ‘stability’ in the sense of the persistence of species [47].
We are now in a position to examine different aspects of stability, including
robustness to secondary extinctions, based on structures directly inferred
from empirical networks. Although the Bayesian modeling approach is not new to
network analysis in general [16, 17], it remains relatively rare. The Bayesian
group model, and, more importantly, the general framework for modeling and
model selection, naturally extend to other kinds of biological networks, such
as metabolic and regulatory networks [52] and networks describing other
ecological interactions such as pollination [53]. We advocate this framework
as a way to build stronger ties between hypothesis formulation, model
building, and data analysis.

## Methods

### Group Model

We use as a starting point the group model of Allesina and Pascual [8], in
which a network of $N$ nodes is partitioned into $K$ groups. The groups to
which a potential prey and to which a potential predator belong completely
determine the probability that a feeding relationship exists between them. The
assignment of species to groups is given by the vector
$\mathbf{G}=(g_{1},\dots,g_{n})$, with $g_{i}\in\\{1,\dots,K\\}$. We refer to
this assignment as a set ‘partition,’ in keeping with standard mathematical
terminology. The probability that a species assigned to group $i$ is consumed
by a species assigned to group $j$ is equal to $p_{ij}$. This gives a matrix
$\mathbf{P}$ of $K^{2}$ probabilities, containing the probabilities of
observing directed links between members of each pair of groups, and within
members of each group.

If we take $\mathbf{A}$ to be the directed adjacency matrix of a network, with
entries $a_{ij}$ equal to 1 if a link exists from node $i$ to node $j$, 0
otherwise, then the probability of the network being generated by partition
$\mathbf{G}$ and link probabilities $\mathbf{P}$ is given by

$\displaystyle f(\mathbf{A}|\mathbf{G},\mathbf{P})$
$\displaystyle=\prod_{i=1}^{K}\prod_{j=1}^{K}p_{ij}^{Y_{ij}}(1-p{ij})^{Z_{ij}}\,,$
(6)

where $Y_{ij}$ and $Z_{ij}$ are the number of 1-entries and 0-entries in the
submatrix of $\mathbf{A}$ containing entries from rows $r$ satisfying
$g_{r}=i$ and columns $c$ satisfying $g_{c}=j$.

In the simplest case, all nodes are assigned to the same group, and the
likelihood simplifies to

$\displaystyle f(\mathbf{A}|p)$ $\displaystyle=p^{Y}(1-p)^{Z}$ (7)

where $Y$ and $Z$ are the total number of 1-entries and 0-entries in
$\mathbf{A}$.

In order to use the group model for Bayesian inference, we want to infer the
posterior distribution over partitions and parameters,

$\displaystyle f(\mathbf{G},\mathbf{P}|\mathbf{A})\propto
f(\mathbf{G},\mathbf{P})f(\mathbf{A}|\mathbf{G},\mathbf{P}).$ (8)

This requires specifying a prior distribution over partitions $\mathbf{G}$ and
link probabilities $\mathbf{P}$. We consider two priors over $\mathbf{G}$ and
two priors over $\mathbf{P}$.

### Model Priors

#### Priors for Partitions

The simplest prior over partitions assigns equal probability to each possible
assignment of nodes into groups. For a network of $N$ nodes, the number of
possible partitions is given by the $N$th Bell number,

$\displaystyle\mathcal{B}(N)$
$\displaystyle=\sum_{K=1}^{N}\mathcal{S}_{2}(N,K),$ (9)

where $\mathcal{S}_{2}(N,K)$ is the Stirling number of the second kind, the
number of ways to partition $N$ objects into exactly $K$ groups,

$\displaystyle\mathcal{S}_{2}(N,K)$
$\displaystyle=\frac{1}{K!}\sum_{j=0}^{K}(-1)^{K-j}{\binom{K}{j}}j^{N}.$ (10)

Therefore, the prior probability of a particular partition is uniform across
all possible partitions

$\displaystyle f(\mathbf{G})$ $\displaystyle=\frac{1}{\mathcal{B}(N)},$ (11)

and the prior probability of having exactly $K$ groups is

$\displaystyle f(K)$
$\displaystyle=\frac{\mathcal{S}_{2}(N,K)}{\mathcal{B}(N)}.$ (12)

For partitions, the choice of a uniform prior, although simple, includes
hidden assumptions. In particular, there are far more possible partitions for
an intermediate number of groups than a small or large number, so the prior
will implicitly bias results toward that number. For example, with 100 nodes,
the distribution is peaked at $K=28$ (Figure S1).

An alternate prior for partitioning objects into groups comes from the
Dirichlet process, also known as the “Chinese restaurant process,” which is
becoming a standard Bayesian prior for related problems [28, 54, 55]. Consider
a restaurant with an infinitely large number of infinitely large tables, all
initially empty. The first patron sits alone at the first table, and
subsequent patrons may either sit at an occupied table or a new table. They
choose occupied tables with weight equal to the number of current occupants,
or a new table with weight equal to an aggregation parameter $\chi$. For
example, the second patron will sit at the same table as the first patron with
probability $1/(1+\chi)$. In fact, because the process is exchangeable, the
probability of any pair of patrons sitting at the table is also $1/(1+\chi)$.
If $\chi$ is small, there will tend to be a small number of occupied tables
and a skewed distribution of table sizes; if $\chi$ is large, there will be a
larger number of tables occupied by few patrons.

Interpreting tables of patrons as groups of nodes, under the Dirichlet process
the prior probability of a particular partition $\mathbf{G}$ is

$\displaystyle f(\mathbf{G}|\chi)$
$\displaystyle=\chi^{K}\frac{\prod_{j=1}^{K}(\eta_{j}-1)!}{\prod_{i=1}^{N}(\chi+i-1)},$
(13)

where $N$ is the number of nodes in the network, $K$ is the number of groups
in the partition, and $\eta_{j}$ is the number of nodes in group $j$. The
prior probability of $K$ groups is

$\displaystyle f(K|\chi)$
$\displaystyle=\frac{|\mathcal{S}_{1}(N,K)|\chi^{K}}{\prod_{i=1}^{N}(\chi+i-1)},$
(14)

where $\mathcal{S}_{1}(N,K)$ is a Stirling number of the first kind, equal to
the coefficients on $x_{K}$ in the expansion $x(x-1)(x-2)\ldots(x-K+1)$.

Rather than choosing a fixed value of $\chi$ for the prior, we give $\chi$ an
exponential hyperprior distribution with mean $1$:

$\displaystyle f(\chi)$ $\displaystyle=e^{-\chi}\qquad\chi\geq 0.$ (15)

#### Priors for Link Probabilities

Similarly, the elements of link probability matrix $\mathbf{P}$ may be given a
simple uniform prior over $[0,1]$:

$\displaystyle f(p_{ij})$ $\displaystyle=1\qquad 0\leq p_{ij}\leq 1.$ (16)

As there may be some regularity in the values of the link probabilities, we
also tried a beta prior:

$\displaystyle f(p_{ij}|\alpha,\beta)$
$\displaystyle=\frac{1}{B(\alpha,\beta)}p_{ij}^{\alpha-1}(1-p_{ij})^{\beta-1},$
(17)

where $B(\alpha,\beta)$ is the beta function,

$\displaystyle B(\alpha,\beta)$
$\displaystyle=\int_{0}^{1}t^{\alpha-1}(1-t)^{\beta-1}\,dt.$ (18)

The parameters $\alpha$ and $\beta$ control the shape of the distribution,
which may be convex, concave, or skewed toward $0$ or $1$. When
$\alpha=\beta=1$, the beta prior becomes a uniform distribution.

We use $\alpha$ and $\beta$ exponential hyperpriors with mean $1$:

$\displaystyle f(\alpha)$ $\displaystyle=e^{-\alpha}\qquad\alpha\geq 0,$ (19)
$\displaystyle f(\beta)$ $\displaystyle=e^{-\beta}\qquad\beta\geq 0.$ (20)

### Markov-chain Monte Carlo Sampling

For networks of any appreciable size, the number of possible partitions is far
too large to enumerate, so we must use a Markov-chain Monte Carlo (MCMC)
technique to sample from the posterior distribution. Here we describe the
sampling procedure and the details of the Metropolis-Hastings proposal
distribution used.

#### Sampling Procedure

We employ the standard Metropolis-Hastings algorithm to sample from the
posterior distribution over partitions and model hyperparameters [56, 57]. The
general idea of an MCMC method is to set up a sequence of dependent samples
$\theta_{1},\theta_{2},\ldots$ that is guaranteed to converge to a target
distribution, in this case the posterior distribution of our model. Starting
from the current sample, a change is proposed, drawn from a proposal
distribution over possible changes, $q(\theta\rightarrow\theta^{*})$. This
sample is either rejected, in which case the current sample is repeated, or
the proposed sample is accepted as the new sample. The Metropolis-Hastings
acceptance probability,

$\displaystyle r(\theta\rightarrow\theta^{*})$
$\displaystyle=\min\left\\{1,\frac{f(\theta^{*})}{f(\theta)}\frac{q(\theta^{*}\rightarrow\theta)}{q(\theta\rightarrow\theta^{*})}\right\\},$

guarantees that the sequence of samples will converge to the posterior
distribution, $f(\theta|D)\propto f(\theta)f(D|\theta)$, the prior times the
likelihood.

For the group model, the samples $\theta$ consist of hyperparameters for the
model variant—the Dirichlet process prior parameter $\chi$ and the beta prior
parameters $\alpha$ and $\beta$—as well as the group count $K$ and assignment
vector $\mathbf{G}$. The link probabilities $\mathbf{P}$ governing links
between groups are not included, because the likelihood function would not be
compatible between partitions with different values of $K$. One possible
solution to this problem would be to include $\mathbf{P}$ in the sampling
procedure, restrict $K$ to a particular number for a particular run, and then
appropriately weight runs with different values of $K$. Another approach is
reversible-jump MCMC [58], which appropriately handles a mapping between two
different parameter spaces as part of the Metropolis-Hastings proposal ratio.
(We tried a reversible-jump scheme, but chains tended to get stuck at local
maxima.)

Instead of trying to sample values of $\mathbf{P}$, we use the marginal
likelihood of a partition given model hyperparameters—that is, the posterior
distribution, conditional on values of $\alpha$, $\beta$, and $\mathbf{G}$,
integrated over all possible values of $\mathbf{P}$—directly in the
Metropolis-Hastings procedure. This is possible because the marginal
likelihood of a single partition can be calculated analytically.

For a beta prior over link probabilities, the likelihood of $\mathbf{G}$,
$\alpha$ and $\beta$ marginalized over all possible values of $\mathbf{P}$ is

$\displaystyle f(\mathbf{A}|\mathbf{G},\alpha,\beta)$
$\displaystyle=\int_{\mathbf{P}}f(\mathbf{P|\alpha,\beta})f(\mathbf{A}|\mathbf{G},\mathbf{P})\,d\mathbf{P}$
(21)
$\displaystyle=\prod_{i=1}^{K}\prod_{j=1}^{K}\int_{0}^{1}\frac{1}{B(\alpha,\beta)}p_{ij}^{\alpha-1}(1-p_{ij})^{\beta-1}p_{ij}^{Y_{ij}}(1-p_{ij})^{Z_{ij}}\,dp_{ij}$
(22)
$\displaystyle=\prod_{i=1}^{K}\prod_{j=1}^{K}\int_{0}^{1}\frac{1}{B(\alpha,\beta)}p_{ij}^{Y_{ij}+\alpha-1}(1-p_{ij})^{Z_{ij}+\beta-1}\,dp_{ij}$
(23)
$\displaystyle=\prod_{i=1}^{K}\prod_{j=1}^{K}\frac{B(Y_{ij}+\alpha,Z_{ij}+\beta)}{B(\alpha,\beta)}.$
(24)

Similarly, for a uniform prior over link probabilities, the marginal
likelihood of a particular partition is simply

$\displaystyle f(\mathbf{A}|\mathbf{G})$
$\displaystyle=\prod_{i=1}^{K}\prod_{j=1}^{K}B(Y_{ij}+1,Z_{ij}+1).$ (25)

#### Proposal Distribution

For the case of uniform priors on partitions and link probabilities, the
proposal distribution only allows changes to the partition. With the Dirichlet
process prior on partitions and the beta prior on link probabilities,
hyperparameters $\chi$, $\alpha$, and $\beta$ can also be changed.

The proposal distribution is described as follows:

  1. 1.

Each hyperparameter $h$ ($\alpha$, $\beta$, and $\chi$) is chosen for update
with probability $p_{h}$, where $p_{h}$ are tuned to improve convergence. A
proposed new value $h^{\prime}$ is drawn from a uniform distribution between
$\max(0,h-r_{h})$ and $h+r_{h}$, where $r_{h}$ is a proposal radius manually
tuned to improve convergence. (A scale-free proposal could easily be used
instead, and may require less tuning.)

  2. 2.

With probability $\left(1-\sum_{h}p_{h}\right)$, a group-change move is
proposed:

    1. (a)

A node $i$ is chosen uniformly at random as the species to be moved.

    2. (b)

Another node $j\neq i$ is chosen uniformly at random.

    3. (c)

If $i$ and $j$ are in different groups, node $i$ is moved into the group of
node $j$. If $i$ and $j$ are in the same group, node $i$ is moved into a new
group.

### Metropolis-coupled MCMC

Although the Metropolis-Hastings algorithm is guaranteed to converge to the
target distribution at some point, local maxima in the likelihood surface can
cause a chain to become stuck for long periods of time. One approach to
avoiding this problem, known as “Metropolis coupling,” involves running
multiple chains in parallel. One chain, the “cold chain,” explores the target
distribution, while the other chains, “hot chains,” explore low-likelihood
configurations more freely. Periodically, swaps are proposed between chains,
allowing good configurations discovered on hot chains to propagate toward the
cold chain.

Rather than exploring the target distribution $f(\theta|D)\propto
f(\theta)f(D|\theta)$, heated chains explore

$\displaystyle f_{\tau}(\theta|D)\propto
f(\theta)\left[f(D|\theta)\right]^{\tau}\qquad\tau\in[0,1],$ (26)

where $\tau$ is a heating parameter. We use linearly spaced values of $\tau$,
with the hottest chain exploring the prior ($\tau=0$) and the coldest chain
exploring the posterior ($\tau=1$).

Swap moves are standard Metropolis-Hastings proposals, but rather than
considering a change to a single chain, they consider a change to the joint
distribution of two chains. The acceptance probability is thus the ratio of
the joint distribution after and before the move:

$\displaystyle
r\left((\theta_{i},\theta_{j})\rightarrow(\theta_{j},\theta_{i})\right)$
$\displaystyle=\frac{f(\theta_{j})\left[f(D|\theta_{j})\right]^{\tau_{i}}f(\theta_{i})\left[f(D|\theta_{i})\right]^{\tau_{j}}}{f(\theta_{i})\left[f(D|\theta_{i})\right]^{\tau_{i}}f(\theta_{j})\left[f(D|\theta_{j})\right]^{\tau_{j}}}$
(27)
$\displaystyle=\left[\frac{f(D|\theta_{i})}{f(D|\theta_{j})}\right]^{\tau_{j}-\tau_{i}},$
(28)

where $\theta_{i},\theta_{j}$ are the configurations that begin in chains $i$
and $j$, and $\tau_{i},\tau_{j}$ are the heat parameters of the two chains.

The use of multiple heated chains has the side effect of drastically improving
estimates of marginal likelihoods for model selection, as described in the
next section.

### Model Selection via Marginal Likelihood

The marginal likelihood of a model is the likelihood averaged over the prior
distribution. That is, it is the likelihood one would expect by randomly
sampling parameters from the prior distribution:

$\displaystyle f(D|M)$
$\displaystyle=\int_{\theta}f(\theta)f(D|\theta)\,d\theta\,.$ (29)

This value serves as a useful measure of model fit because it directly
incorporates the dependence of the likelihood on uncertainty in parameter
values, implicitly penalizing extra degrees of freedom [49]. If an additional
parameter improves the maximum likelihood but decreases the average
likelihood, the model suffers from overfitting relative to the simpler model.

#### Marginal Likelihood for the Group Model

The marginal likelihood for the group model involves integrating over all
hyperparameters, partitions, and link probabilities. For the model with
uniform distributions over partitions and link probabilities, the marginal
likelihood is

$\displaystyle f(\mathbf{A}|M_{u,u})$
$\displaystyle=\sum_{\mathbf{G}}f(\mathbf{G})f(\mathbf{A}|\mathbf{G})$ (30)
$\displaystyle=\sum_{\mathbf{G}}\frac{1}{\mathcal{B}(N)}\left[\prod_{i=1}^{K}\prod_{j=1}^{K}B(Y_{ij}+1,Z_{ij}+1)\right].$
(31)

With a Dirichlet process prior over partitions and a uniform distribution over
link probabilities, the marginal likelihood is similarly

$\displaystyle f(\mathbf{A}|M_{d,u})$
$\displaystyle=\int_{0}^{\infty}f(\chi)\sum_{\mathbf{G}}f(\mathbf{G}|\chi)f(\mathbf{A}|\mathbf{G})\,d\chi.$
(32)

Using a uniform prior over partitions and a beta prior over link probabilities
yields

$\displaystyle f(\mathbf{A}|M_{u,b})$
$\displaystyle=\sum_{\mathbf{G}}f(\mathbf{G})\int_{0}^{\infty}f(\alpha)\int_{0}^{\infty}f(\beta)f(\mathbf{A}|\mathbf{G},\alpha,\beta)\,d\beta\,d\alpha.$
(33)

Combining both gives

$\displaystyle f(\mathbf{A}|M_{d,b})$
$\displaystyle=\int_{0}^{\infty}f(\chi)\sum_{\mathbf{G}}f(\mathbf{G}|\chi)\int_{0}^{\infty}f(\alpha)\int_{0}^{\infty}f(\beta)f(\mathbf{A}|\mathbf{G},\alpha,\beta)\,d\beta\,d\alpha\,d\chi.$
(34)

#### Thermodynamic Integration for Marginal Likelihood Estimation

As enumeration across all possible partitions is impossible for networks of
any significant size, we would like to use MCMC to estimate the marginal
likelihood for the sake of comparison among different models. Marginal
likelihood estimates derived from a single chain, such as the harmonic mean
estimator of Raftery [31], converge very slowly, because MCMC fails to sample
sufficiently from low-likelihood areas. However, it is possible to use the
information gathered about low-likelihood areas in heated chains using a
technique called thermodynamic integration [59, 60], or path sampling [61].

Assuming a continuum of heated chains, the thermodynamic estimator for the
log-marginal likelihood is

$\displaystyle\log\hat{\mathcal{L}}(M)$
$\displaystyle=\int_{0}^{1}\frac{1}{m}\sum_{i=1}^{m}\pi(\theta_{i,\tau})\log\mathcal{L}(\theta_{i,\tau})\,d\tau$
(35)

where $m$ is the number of samples in the MCMC output, and $\theta_{i,\tau}$
is a single sample from the output in a chain with heat parameter $\tau$ [60].
With a finite number of chains, we estimate this integral using cubic spline
interpolation as implemented in the splinefun function in the R software
package [62].

### Consensus Partitions

The full output of an MCMC chain from the group model includes an extremely
large number of different partitions, and, for the sake of interpretation, it
is desirable to seek a consensus partition that does a reasonable job of
summarizing the distribution. We use a simple, computationally inexpensive
method to accomplish this task: in short, clustering the nodes in the network
based on a pairwise group-membership matrix.

The group-membership matrix $\mathbf{M}$ is the posterior probability that two
nodes are in the same group and 0 otherwise, that is,

$\displaystyle\mathbf{M}$
$\displaystyle=\sum_{\mathbf{G}}P(\mathbf{G}|\mathbf{A})\mathbf{M}_{\mathbf{G}}\,,$
(36)

where an entry $\mathbf{M}_{\mathbf{G}}$ is 1 if nodes $i$ and $j$ are in the
same group, that is,

$\displaystyle\mathbf{M}_{\mathbf{G},ij}$
$\displaystyle=\delta_{\mathbf{G}_{i},\mathbf{G}_{j}}\,,$ (37)

where $\delta$ is the Kronecker delta and $\mathbf{G}$ is the assignment
vector for the partition. This matrix is estimated from MCMC output as the
fraction of MCMC samples in which the corresponding species are in the same
group:

$\displaystyle\mathbf{\hat{M}}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\mathbf{M}_{\mathbf{G}_{i}}\,.$ (38)

A consensus partition is formed by applying a hierarchical clustering
algorithm to the group-membership matrix estimate $\mathbf{\hat{M}}$, and then
cutting the dendrogram at some number of groups $K$, forming a consensus
partition with assignment vector $\mathbf{G}_{K}$ and group-membership matrix
$\mathbf{M}_{K}$. The goodness of fit of a consensus partition is simply
measured as the correlation between $\mathbf{\hat{M}}$ and $\mathbf{M}_{K})$.
The best consensus partition is thus identified using the value of $K$ that
gives the highest correlation.

We use the average-linkage clustering algorithm [63] as implemented by the
hclust function in the R software package [62], treating
$\mathbf{1}-\mathbf{\hat{M}}$ as distance matrix. We find that the average-
linkage algorithm produces higher correlations than the other algorithms
implemented as well as ideal $K$ close to the mean $K$ in the MCMC output.
Furthermore, we find that consensus partitions produce higher correlations
with the $\mathbf{\hat{M}}$ than any individual partition in the MCMC output.

## Acknowledgments

We acknowledge the support of NSF grant EF-0827493 (Program on Theory in
Biology) to S.A. and M.P., and of the DOE Computational Science Graduate
Fellowship (grant DE-FG02-97ER25308) to E.B. T.B. was supported by the Howard
Hughes Medical Institute, and M.P. is a Howard Hughes Medical Institute
Investigator.

A.D. acknowledges the McDonnell Foundation for financial support, the
Frankfurt Zoological Society for logistical support in the Serengeti for his
work on food webs, and the members of Serengeti Biocomplexity Project for many
interesting discussions about the Serengeti food web.

The funders had no role in study design, data collection and analysis,
decision to publish, or preparation of the manuscript.

## References

  * 1. Cohen JE (1995) Frontiers in Mathematical Biology, Springer, volume 100 of _Lecture Notes in Biomathematics_.  pp. 351–380.
  * 2. Camerano L (1880) Dell’equilibrio dei viventi merce la reciproca distruzione.  Accademia delle Scienze di Torino 15: 393–414.
  * 3. Elton CS (1927) Animal Ecology.  Text-Books of Animal Biology. London: Sidgwick & Jackson, Ltd.
  * 4. May RM (1973) Stability and Complexity in Model Ecosystems.  Princeton University Press.
  * 5. Pimm S, Lawton J (1980) Are food webs divided into compartments?  The Journal of Animal Ecology 49: 879–898.
  * 6. Krause A, Frank K, Mason D, Ulanowicz R, Taylor W (2003) Compartments revealed in food-web structure.  Nature 426: 282–285.
  * 7. Guimera R, Stouffer D, Sales-Pardo M, Leicht E, Newman M, et al. (2010) Origin of compartmentalization in food webs.  Ecology 91: 2941–2951.
  * 8. Allesina S, Pascual M (2009) Food web models: a plea for groups.  Ecology Letters 12: 652–62.
  * 9. Girvan M, Newman M (2002) Community structure in social and biological networks.  Proceedings of the National Academy of Sciences 99: 7821.
  * 10. Burns T (1989) Lindeman’s contradiction and the trophic structure of ecosystems.  Ecology 70: 1355–1362.
  * 11. Luczkovich JJ, Borgatti SP, Johnson JC, Everett MG (2003) Defining and measuring trophic role similarity in food webs using regular equivalence.  Journal of Theoretical Biology 220: 303–321.
  * 12. Cohen JE, Briand F, Newman CM (1990) Community Food Webs: Data and Theory.  Biomathematics. Springer.
  * 13. Williams R, Martinez N (2000) Simple rules yield complex food webs.  Nature 404: 180–183.
  * 14. Allesina S, Alonso D, Pascual M (2008) A general model for food web structure.  Science 320: 658–61.
  * 15. McCarthy MA (2007) Bayesian Methods for Ecology.  Cambridge University Press.
  * 16. Hoff P, Raftery A, Handcock M (2002) Latent space approaches to social network analysis.  Journal of the American Statistical Association 97: 1090–1098.
  * 17. Park Y, Moore C, Bader JS (2010) Dynamic networks from hierarchical Bayesian graph clustering.  PLoS One 5: e8118.
  * 18. Clauset A, Moore C, Newman MEJ (2008) Hierarchical structure and the prediction of missing links in networks.  Nature 453: 98–101.
  * 19. Williams R, Anandanadesan A, Purves D (2010) The probabilistic niche model reveals the niche structure and role of body size in a complex food web.  PLoS One 5: e12092.
  * 20. Yang Z, Rannala B (1997) Bayesian phylogenetic inference using DNA sequences: a Markov chain Monte Carlo method.  Molecular Biology and Evolution 14: 717.
  * 21. Mau B, Newton M, Larget B (1999) Bayesian phylogenetic inference via Markov chain Monte Carlo methods.  Biometrics 55: 1–12.
  * 22. Dobson A (2009) Food-web structure and ecosystem services: insights from the Serengeti.  Philosophical Transactions of the Royal Society B: Biological Sciences 364: 1665–1682.
  * 23. Sinclair ARE, Norton-Griffiths M, editors (1979) Serengeti, Dynamics of an Ecosystem.  University Of Chicago Press.
  * 24. Sinclair ARE, Arcese P, editors (1995) Serengeti II: Dynamics, Management, and Conservation of an Ecosystem.  University Of Chicago Press.
  * 25. Sinclair ARE, Packer C, Mduma SA, Fryxell JM, editors (2008) Serengeti III: Human Impacts on Ecosystem Dynamics.  University Of Chicago Press.
  * 26. McNaughton S (1983) Serengeti grassland ecology: the role of composite environmental factors and contingency in community organization.  Ecological Monographs 53: 291–320.
  * 27. Wang Y, Wong G (1987) Stochastic blockmodels for directed graphs.  Journal of the American Statistical Association 82: 8–19.
  * 28. Ferguson T (1973) A Bayesian analysis of some nonparametric problems.  The Annals of Statistics 1: 209–230.
  * 29. Jeffreys H (1935) Some tests of significance, treated by the theory of probability.  Proceedings of the Cambridge Philosophical Society 31: 203–222.
  * 30. Jeffreys H (1961) Theory of Probability.  The International Series of Monographs on Physics. Oxford: Clarendon Press, 3rd edition.
  * 31. Kass R, Raftery A (1995) Bayes factors.  Journal of the American Statistical Association 90.
  * 32. Burnham KP, Anderson D (2002) Model Selection and Multi-Model Inference.  Springer.
  * 33. Sinclair A, Mduma S, Brashares J (2003) Patterns of predation in a diverse predator–prey system.  Nature 425: 288–290.
  * 34. Casebeer R, Koss G (1970) Food habits of wildebeest, zebra, hartebeest and cattle in Kenya Masailand.  African Journal of Ecology 8: 25–36.
  * 35. Cooper S, Holekamp K, Smale L (1999) A seasonal feast: long-term analysis of feeding behaviour in the spotted hyaena (Crocuta crocuta).  African Journal of Ecology 37: 149–160.
  * 36. Hansen R, Mugambi M, Bauni S (1985) Diets and trophic ranking of ungulates of the northern Serengeti.  The Journal of Wildlife Management : 823–829.
  * 37. Murray M (1993) Comparative nutrition of wildebeest, hartebeest and topi in the Serengeti.  African Journal of Ecology 31: 172–177.
  * 38. Talbot L, Talbot M (1962) Food preferences of some East African wild ungulates.  East African Agricultural and Forestry Journal 27: 131–138.
  * 39. Talbot LM, Talbot MH (1963) The Wildebeest in Western Masailand, East Africa.  Wildlife Monographs : 3–88.
  * 40. Schaller G (1972) The Serengeti Lion.  University of Chicago Press.
  * 41. Caro T (1994) Cheetahs of the Serengeti Plains: Group Living in an Asocial Species.  University of Chicago Press.
  * 42. Kruuk H (1972) The Spotted Hyaena.  University of Chicago Press.
  * 43. McNaughton S (1985) Ecology of a grazing ecosystem: the Serengeti.  Ecological Monographs 55: 259–294.
  * 44. Vesey-FitzGerald D (1960) Grazing succession among East African game animals.  Journal of Mammalogy 41: 161–172.
  * 45. Lamprecht J (1978) On diet, foraging behaviour and interspecific food competition of jackals in the Serengeti National Park, East Africa.  Zeitschrift für Säugetierkunde 43: 210–223.
  * 46. Rooney N, McCann K, Moore J (2008) A landscape theory for food web architecture.  Ecology Letters 11: 867–881.
  * 47. McCann K, Rasmussen J, Umbanhowar J (2005) The dynamics of spatially coupled food webs.  Ecology Letters 8: 513–523.
  * 48. Rezende E, Albert E, Fortuna M, Bascompte J (2009) Compartments in a marine food web associated with phylogeny, body mass, and habitat structure.  Ecology Letters 12: 779–788.
  * 49. Bolker BM (2008) Ecological Models and Data in R.  Princeton University Press.
  * 50. Good B, De Montjoye Y, Clauset A (2010) Performance of modularity maximization in practical contexts.  Physical Review E 81: 46106.
  * 51. Newman M (2006) Modularity and community structure in networks.  Proceedings of the National Academy of Sciences 103: 8577.
  * 52. Jeong H, Tombor B, Albert R, Oltvai Z, Barabási A (2000) The large-scale organization of metabolic networks.  Nature 407: 651–654.
  * 53. Bascompte J, Jordano P, Melián C, Olesen J (2003) The nested assembly of plant–animal mutualistic networks.  Proceedings of the National Academy of Sciences 100: 9383.
  * 54. Huelsenbeck JP, Suchard MA (2007) A nonparametric method for accommodating and testing across-site rate variation.  Systematic Biology 56: 975–87.
  * 55. Xing E, Jordan M, Sharan R (2007) Bayesian haplotype inference via the Dirichlet process.  Journal of Computational Biology 14: 267–284.
  * 56. Metropolis N, Rosenbluth A, Rosenbluth M, Teller A, Teller E (1953) Equation of state calculations by fast computing machines.  Journal of Chemical Physics 21: 1087.
  * 57. Hastings W (1970) Monte Carlo sampling methods using Markov chains and their applications.  Biometrika 57: 97–109.
  * 58. Green P (1995) Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.  Biometrika 82: 711.
  * 59. Lartillot N, Philippe H (2006) Computing Bayes factors using thermodynamic integration.  Systematic Biology 55: 195–207.
  * 60. Beerli P, Palczewski M (2010) Unified framework to evaluate panmixia and migration direction among multiple sampling locations.  Genetics 185: 313–26.
  * 61. Gelman A, Meng X (1998) Simulating normalizing constants: from importance sampling to bridge sampling to path sampling.  Statistical Science 13: 163–185.
  * 62. R Development Core Team (2010) R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.  URL http://www.R-project.org.  ISBN 3-900051-07-0.
  * 63. Sokal R, Michener C (1958) A statistical method for evaluating systematic relationships.  University of Kansas Science Bulletin 38: 1409–1438.
  * 64. Duong T (2007) ks: Kernel density estimation and kernel discriminant analysis for multivariate data in R.  Journal of Statistical Software 21: 1–16.

## Figures

Figure 1: The Serengeti food web. The network is shown using a spring-layout
algorithm without clustering. Plant, herbivore, and carnivore nodes are green,
orange, and red, respectively.  Figure 2: Posterior distributions and prior
expectations of aggregation parameter $\chi$ and group count $K$. Figure 3:
Pairwise group membership matrix. Species are identically ordered top to
bottom and left to right according to the $16$-group consensus partition as
listed in Table 2. Hue indicates group identity; color saturation indicates
the fraction of partitions in which species occupy the same group. Note that
this image conveys information about group membership, not network
connectivity. Figure 4:  Adjacency matrix ordered by groups. Species are
identically ordered top to bottom and left to right according to the
$16$-group consensus partition as listed in Table 2. Black matrix entries
indicate that the species in the column feeds on the species in the row.
Columns that would indicate prey of plant groups are omitted. Note that in a
modular network according to the standard definition, links would be
concentrated on the diagonal of the adjacency matrix, since they occur within
groups. By contrast, here links are concentrated in off-diagonal blocks.
Figure 5:  Network layout of groups. The network is shown organized and
colored by group according to the $16$-group consensus partition listed in
Table 2. Figure 6:  Network layout of aggregated groups. Nodes in the network
are aggregated and colored by group according to the $16$-group consensus
partition listed in Table 2, and arranged vertically by trophic level. Line
thickness indicates the link density between groups. Node area is proportional
to the number of species in a group.

## Supporting Information

Figure S1: Implicit prior on number of groups with uniform partition prior.
The prior distribution on number of groups $K$ is shown for a uniform
partition prior for a network with 100 nodes. The mode of the distribution is
at $K=28$. Figure S2: Posterior distributions of link density parameters
$\alpha$ and $\beta$. Color brightness indicates posterior density, estimated
using the ks multivariate kernel density estimation package for R [64].
Contours indicate cumulative density. The $\alpha$ parameter is significantly
lower than 1, indicating departure from a uniform distribution. Figure S3:
Distribution of link probability parameters. The prior distribution for link
probability parameters, integrated over the priors for beta distribution
parameters $\alpha$ and $\beta$, is indicated with a dotted line. The heat map
shows beta distributions corresponding to the posterior distribution for
$\alpha$ and $\beta$, with lightness indicating the posterior density of the
parameter values.

(See supporting file tab_S1_species_list.csv for data.)

Table S1: Species in the Serengeti food web.

(See supporting file tab_S2_edge_list.csv for data.)

Table S2: Feeding links in the Serengeti food web.

(See supporting file tab_S3_consensus_partition.csv for data.)

Table S3: $\mathbf{16}$-group consensus partition.

(See supporting file tab_S4_link_density.csv for data.)

Table S4: Link densities between groups in the $\mathbf{16}$-group consensus
partition.

## Tables

Table 1:  Marginal likelihood estimates for model variants calculated via
thermodynamic integration.

Partition prior | Link prior | Log marginal likelihood estimate | 95% bootstrap confidence int.
---|---|---|---
Uniform | Uniform | $-1826.87$ | $(-1826.94,-1826.82)$
Uniform | Beta | $-1463.22$ | $(-1463.39,-1463.03)$
Dirichlet process | Uniform | $-1547.41$ | $(-1547.45,-1547.37)$
Dirichlet process | Beta | $-1356.96$ | $(-1357.02,-1356.90)$
One group | Uniform | $-2870.58$ | (exact)
$170$ groups | Uniform | $-20031.95$ | (exact)
$170$ groups | Beta | $-2870.58$ |

Table 2:  Groups identified in the Serengeti food web using a $\mathbf{16}$-group consensus partition. Group 1 |  Crocuta crocuta, Lycaon pictus, Panthera leo, Panthera pardus, Acinonyx jubatus
---|---
Group 2 |  Canis aureus
Group 3 |  Canis mesomelas, Leptailurus serval, Caracal caracal
Group 4 |  Connochaetus taurinus, Gazella granti, Gazella thomsoni, Equus burchelli, Alcelaphus buselaphus, Aepyceros melampus, Damaliscus korrigum
Group 5 |  Kobus ellipsiprymnus, Phacochaerus aethiopicus, Tragelaphus scriptus, Ourebia ourebi, Redunca redunca, Pedetes capensis, Taurotragus oryx, Rhabdomys pumilio, Hippopotamus amphibus, Cercopithecus aethiops
Group 6 |  Syncerus caffer
Group 7 |  Heterohyrax brucei, Procavia capensis
Group 8 |  Agama planiceps
Group 9 |  Papio anubis, Giraffa camelopardalis, Madoqua kirkii
Group 10 |  Loxodenta africana
Group 11 |  Panicum coloratum, Sporobolus pyramidalis, Hyparrhenia filipendula, Harpachne schimperi, Digitaria macroblephara, Eragrostis tenuifolia, Grewia bicolor, Aristida adoensis, Brachiaria semiundulata, Pennisetum mezianum, Bothriochloa insculpta, Panicum maximum, Sida spp., Eustachys paspaloides, Croton macrostachyus, Solanum incanum, Indigofera hochstetteri, Hibiscus spp., Heteropogon contortus, Cynodon dactylon, Themeda triandra, Balanites aegytiaca
Group 12 |  Digitaria scalarum/abysinnica, Dinebra retroflexa, Ischaemum afrum, Eragostris cilianensis, Hyparrhenia rufa, Sporobolus fimbriatus, Sporobolus spicatus
Group 13 |  Microchloa kunthii
Group 14 |  Echinochloa haploclada, Digitaria milanjiana, Panicum deustum, Digitaria ternata, Andropogon schirensis, Cymbopogon excavatus, Setaria sphacelata, Typha capensis, Setaria pallidifusca, Phragmites mauritianus, Eragrostis exasperata, Andropogon greenwayi, Lonchocarpus eriocalyx, Sporobolus centrifugus, Hyparrhenia dissoluta, Chloris roxburghiana, Aristidia hordacea, Chloris pycnothrix, Panicum repens, Aristidia kenyensis, Combretum molle, Acacia xanthophloea, Disperma kilimandscharica, Vossia cuspida, Odyssea jaegeri, Sporobolus ioclados, Euphorbia candelabrum, Sorghum versicolor, Kigelia africana, Olea spp., Sporobolus festivus, Acacia pallens, Crotalaria spinosa, Digitaria diagonalis, Boscia augustifolia, Acacia robusta, Acacia seyal/hockii, Chloris gayana, Pennisetum stramineum, Commiphora africana/trothae
Group 15 |  Acacia senegal, Acacia tortilis
Group 16 |  Grewia fallax, Cissus quadrangularis, Cissus rotundifolia, Commelina africana, Allophylus rubifolus, Sensevieria ehrenbergiana, Pavetta assimilis, Phyllanthus sepialis, Acalypha fructicosa, Maerua triphllya, Ficus glumosa, Croton dichogamus, Sclerocarya birrea, Capparis tomentosa, Ximenia caffra, Cordia ovalis, Grewia trichocarpa, Abutilon angulatum, Pappaea capensis, Commiphora schimperi, Albuca spp., Ficus ingens, Hoslundia opposita, Ocinum suave, Cenchrus ciliaris, Solanum dennekense, Aloe macrosiphon, Indigofera basiflora, Ipomoea obscura, Albizia harveyi, Ficus thinningii, Emilia coccinea, Cyphostemma nierensis, Spirocarpa spp., Sensevieria suffruticosa, Pupalia lappacea, Aloe secundiflora, Turreae fischeri, Pavonia patens, Jasminum fluminense, Acacia clavigera, Cassia didymobotrya, Kedrotis foetidissima, Hypoestes forskalii, Zisiphus mucronata, Commiphora merkeri, Blepharis acanthoides, Iboza sp., Rhoicissus revoilii, Kalanchoe sp., Solanum nigrum, Achyranthes aspera, Digitaria velutina, Tricholaena eichingeri, Lippia ukambensis, Heliotropium steudneri, Kyllinga nervosa, Sporobolus stapfianus, Cyperus kilimandscharica, Pellaea calomelanos, Sporobolus pellucidus, Eragrostis aspera, Eriochloa nubica, Diheteropogon amplectus