Saint- Clair Chabert-Liddell

Pierre Barbillon

Sophie Donnet

Learning common structures in a collection of networks. An application to food webs

Let a collection of networks represent interactions within several (social or ecological) systems. We pursue two objectives: identifying similarities in the topological structures that are held in common between the networks and clustering the collection into sub-collections of structurally homogeneous networks. We tackle these two questions with a probabilistic model based approach. We propose an extension of the Stochastic Block Model (SBM) adapted to the joint modeling of a collection of networks. The networks in the collection are assumed to be independent realizations of SBMs. The common connectivity structure is imposed through the equality of some parameters.

The model parameters are estimated with a variational Expectation-Maximization (EM) algorithm. We derive an ad-hoc penalized likelihood criterion to select the number of blocks and to assess the adequacy of the consensus found between the structures of the different networks. This same criterion can also be used to cluster networks on the basis of their connectivity structure. It thus provides a partition of the collection into subsets of structurally homogeneous networks.

The relevance of our proposition is assessed on two collections of ecological networks. First, an application to three stream food webs reveals the homogeneity of their structures and the correspondence between groups of species in different ecosystems playing equivalent ecological roles. Moreover, the joint analysis allows a finer analysis of the structure of smaller networks. Second, we cluster 67 food webs according to their connectivity structures and demonstrate that five mesoscale structures are sufficient to describe this collection.

Introduction

Context The last few years have seen an increase in the number of interaction networks collected, as networks are popular tools for representing the functioning of a social or ecological system. For a long time the statistical analysis of network data has focused on analyzing a single network at a time. This can be performed either by looking at local or global topological features, or by setting a probabilistic model inferred from the network [START_REF] Kolaczyk | Statistical Analysis of Network Data: Methods and Models[END_REF]. When several networks describing the same kind of interactions are available, a natural question is to assess to what extent they are similar or different. As network data are complex by nature, this comparison of different networks is not an easy task and has mainly focused on comparing statistical topological features on the local, global or mesoscale levels. These comparison metrics depend on whether the networks are defined on the same set of nodes or on different sets of nodes . A survey for the former case is dealt with in [START_REF] Donnat | Tracking network dynamics: A survey using graph distances[END_REF] and another survey for both cases is done in [START_REF] Wills | Metrics for graph comparison: a practitioner's guide[END_REF].

In this paper, we consider networks with no node correspondence and no link between networks as it may be the case in multilayer networks [START_REF] Kivelä | Multilayer networks[END_REF]. Furthermore, the networks are assumed to represent interactions of the same type (directed or not) and with the same valuation (binary, discrete or continuous). A set of such networks constitutes what we call in this paper a collection of networks, although some authors may use this terminology in a different meaning.

When observing such a collection, we aim to determine if the respective structures of the networks are similar. This paper focuses on the mesoscale structure of the networks by assuming that the nodes can be grouped into blocks on the basis of their connectivity pattern [START_REF] White | Social structure from multiple networks. i. blockmodels of roles and positions[END_REF]. A classical tool to infer such a mesoscale structure of a single network is the Stochastic Block Model [START_REF] Holland | Stochastic blockmodels: First steps[END_REF]Snijders and Nowicki, 1997, SBM). In the SBM, a latent variable is associated with each node giving its group/block membership. Nodes belonging to the same block share the same connectivity pattern. The SBM has easily interpretable parameters and its framework allows multiple extensions such as modeling the interactions with various distributions [START_REF] Mariadassou | Uncovering latent structure in valued graphs: a variational approach[END_REF]. The block memberships are not known a priori, they are recovered a posteriori by the inference algorithm.

In social (resp. ecological) networks, individuals (resp. species) with the same block membership play the same social/ecological role in its system [START_REF] Boorman | Social Structure from Multiple Networks. II. Role Structures[END_REF]Luczkovich et al., 2003). In food webs, species playing the same ecological role are said to be ecologically equivalent (see [START_REF] Cirtwill | A review of species role concepts in food webs[END_REF], for a review of species role concepts in food webs). When analysing the roles in food webs, Luczkovich et al. (2003) use the notion of regular equivalence to define trophic role. Two species are said to be regularly equivalent if they feed on equivalent species and are preyed on by equivalent species. This notion of regular equivalence is a relaxation of structural equivalence which imposes that structurally equivalent species have exactly the same trophic relations in the food web. In practice, Luczkovich et al. (2003) find that species are grouped into blocks by trophic level and some separation might occur based on trophic chains. Other papers lead to a similar interpretation of the blocks for stochastic equivalence when fitting SBMs on food webs. It is also noticed that communities (blocks of species preying on each other) are unusual [START_REF] Allesina | Food web models: a plea for groups[END_REF][START_REF] Sander | What can interaction webs tell us about species roles?[END_REF]. Stochastic equivalence has the advantage of taking into account the noisy aspects of the observed networks. In addition, SBM, as a probabilistic generative model, provides the modeler with a unified framework for model selection, link prediction, simulation, and modeling extension, for example, for a collection of networks.

Inferring independently an SBM for each network and comparing them may be misleading. Indeed, a given network may have several possible grouping of the nodes into blocks [START_REF] Peel | The ground truth about metadata and community detection in networks[END_REF] that are equally likely [START_REF] Peixoto | Hierarchical block structures and high-resolution model selection in large networks[END_REF]. Furthermore, the observation of a network may be noisy [START_REF] Guimerà | Missing and spurious interactions and the reconstruction of complex networks[END_REF], especially for ecological networks, the sampling of which is known to be incomplete [START_REF] Rivera-Hutinel | Effects of sampling completeness on the structure of plant-pollinator networks[END_REF].

Our contribution Thus, we propose to jointly model a collection of networks by extending the SBM. We assume that the networks are independent realizations of SBMs sharing common parameters. The natural and interesting consequence is the correspondence between the blocks of the different networks. The proposed model called colSBM comes with a few variants. The simplest model assumes that the parameters of the SBMs are identical leading to a collection of i.i.d. networks. As this assumption might be too restrictive for real networks, we introduce two relaxations on this assumption. The first one is to allow the distribution of the block memberships to vary between networks and even to allow some networks to not populate certain blocks. This enables to model a collection of networks where the structure of certain networks is encompassed in the structure of other networks. The second relaxation allows networks to have the same structure up to a density parameter. This is particularly useful to model networks with different sampling efforts, since it has a direct impact on the density of ecological networks [START_REF] Blüthgen | Measuring specialization in species interaction networks[END_REF].

The inference of the block memberships, the model parameters and the model selection are done through an ad-hoc version of classic tools when inferring SBM, namely a Variational EM algorithm for the inference and an adaptation of the integrated classification likelihood (ICL) criterion for the model selection [START_REF] Daudin | A mixture model for random graphs[END_REF].

The interest of our colSBM model is two-folds. The first one is to find a common connectivity pattern which explains the structure of the different networks in the collection and to assess via model selection whether these structures are a reasonable fit for the collection. As a by-product, it allows a fine analysis of the role structure of the different nodes in the networks. By sharing the blocks between the networks, colSBM allows to recover sets of nodes which play the same sociological/ecological roles in different networks. The second one is to provide a partition of the collection of networks into sub-collections of structurally homogeneous networks1 . Both aspects have some practical implications in ecology [START_REF] Ohlsson | Spatial resolution and location impact group structure in a marine food web[END_REF][START_REF] Michalska-Smith | Telling ecological networks apart by their structure: A computational challenge[END_REF].

As a side effect, by modeling these networks together, provided that the networks have common connectivity patterns, we can use the information of certain networks to recover noisy information from other networks by improving the prediction of missing links [START_REF] Clauset | Hierarchical structure and the prediction of missing links in networks[END_REF]. Hence colSBM has a stabilizing effects on the grouping of the nodes into blocks and might give a block membership that is closer to the one of the full real network than just a single SBM as this will be shown in the numerical studies and application.

Related work Since the SBM is a very flexible model, it has already been adapted to multilayer networks. To name a few, [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF] model a collection of networks along a time gradient, the connectivity structure varies from time to time but they integrate a sparsity parameter, which is similar to our density parameter in the binary case.When dealing with networks with no common nodes, [START_REF] Chabert-Liddell | A stochastic block model approach for the analysis of multilevel networks: An application to the sociology of organizations[END_REF] deal with multilevel networks where the networks are linked by a hierarchical relation between the nodes of the different levels. Within the SBM framework, the closest work to ours is the strata multilayer SBM [START_REF] Stanley | Clustering network layers with the strata multilayer stochastic block model[END_REF], in that it looks for both common connectivity patterns and network clustering. However, it does not consider a collection of networks but a multiplex network where all the networks share the same nodes.

Most contributions about collections of networks rely on some node correspondence between the networks. Recently, motivated by the analysis of fMRI data a few works extend the SBM to model population of networks [START_REF] Paul | A random effects stochastic block model for joint community detection in multiple networks with applications to neuroimaging[END_REF][START_REF] Pavlović | Multisubject Stochastic Blockmodels for adaptive analysis of individual differences in human brain network cluster structure[END_REF]. [START_REF] Le | Estimating a network from multiple noisy realizations[END_REF] make the assumption that the networks of the collection are noisy realizations of the true network, while [START_REF] Reyes | Stochastic blockmodels for exchangeable collections of networks[END_REF] use in a Bayesian framework a hierarchical SBM to model the collection. [START_REF] Durante | Nonparametric bayes modeling of populations of networks[END_REF] propose a mixture of latent space models and Signorelli and Wit (2020) a mixture of network models which is not restricted to the SBM.

Dealing with networks with no node correspondence, Faust and Skvoretz (2002) compare networks involving different species and interaction types using the parameters of exponential random graph models (ERGMs or p * models). More recently, [START_REF] Yin | Finite Mixtures of ERGMs for Modeling Ensembles of Networks[END_REF] propose a mixture of ERGMs to model the generative process of a collection of networks. ERGMs allow testing the significance of selected local interaction patterns that convey ecological or sociological meaning. Compared to colSBM , they do not group nodes into blocks and as such do not provide role equivalence between nodes of different networks. The contributions dealing with networks with no node correspondence also include a hierarchical mixed membership SBM, using a common Bayesian prior on the connectivity parameter of the different networks [START_REF] Sweet | Hierarchical mixed membership stochastic blockmodels for multiple networks and experimental interventions[END_REF]. Finally on partitioning a collection of networks, [START_REF] Mukherjee | On clustering network-valued data[END_REF] use graph moments (they also propose to fit a mixture of graphon when having access to node correspondence between the networks in the collection and then to make a spectral clustering on the distance matrix between networks), while [START_REF] Sweet | Clustering ensembles of social networks[END_REF] use graph kernel methods on networks with nodes label to estimate independently a feature vector for each network. Both of these contributions rely on clustering those feature vectors, and as such do not provide any estimate on the joint structure of the collection.

Outline Section 2 recalls the definition of the Stochastic Block Model on a single network. We motivate our new approach by inferring it independently on a collection of food webs. Then in Section 3, we present the various variants of the colSBM . The likelihood expression is provided in Section 4, together with some identifiability conditions. We develop the methodology for the parameter estimation and model selection in Section S-2, while Section 6 deals with network clustering. The details on the clustering procedure are postponed to Appendix A.. We finally propose with two applications on food webs in Section 7. First, we compare the structures of 3 networks and show the information transfer between these networks. Second, we seek a partition of a collection of 67 networks. The technical details and numerical studies which demonstrates the efficiency of our inference procedure and the pertinence of our model selection criterion are left in Supplementary Material [START_REF] Chabert-Liddell | Supplement to "learning common structures in a collection of networks[END_REF].

Data motivation and the stochastic block model

Consider a collection of M independent networks where each network indexed by m involves its own n m nodes. The networks are encoded into their adjacency matrices

(X m ) m∈{1,...,M } such that: ∀m ∈ {1, . . . , M }, ∀(i ̸ = j) ∈ {1, . . . , n m } 2 , X m ij = 0 if no interaction is observed between species i and j of network m X m ij ̸ = 0 otherwise.
If the networks represent binary interactions then X m ij ∈ K = {0, 1}, ∀(m, i, j); if the interactions are weighted such as counts, then X m ij ∈ K = N. Moreover, all the networks encompass the same type of interactions (binary, count. . . ) and no self-interaction is considered. Besides, for the sake of simplicity, we assume that all the networks are directed. The extension to undirected networks i.e. such that X m ij = X m ji for any i ̸ = j is straightforward. X = (X 1 , . . . , X M ) denotes the collection of adjacency matrices.

A first ecological example: three stream food webs As a first example, we consider the collection of three stream food webs from [START_REF] Thompson | Impacts on stream food webs of native and exotic forest: an intercontinental comparison[END_REF]. The three networks collected respectively in Martins (Maine USA), Cooper and Herlzier (North-Carolina, USA) involve respectively 105, 58 and 71 species resulting in 343, 126 and 148 binary edges respectively. Classically, the food web edges represent directed trophic links showing the energy flow ie. X m ij = 1 if species j preys on species i, with no reciprocal interactions. When aiming at unraveling the structure of these networks the Stochastic Block Model (SBM) is an interesting tool which has proven its high flexibility by encompassing a large variety of structures (see [START_REF] Allesina | Food web models: a plea for groups[END_REF], for the particular case of food webs). When dealing with three networks, the standard strategy that we describe below, is to fit separately one SBM per network.

Separate SBM (sepSBM) The SBM introduces blocks of nodes and assumes that the interaction between two nodes is driven by the blocks the nodes belong to. More precisely, for network m, let the n m nodes be divided into Q m blocks. Let Z m = (Z m 1 , . . . , Z m nm ) be independent latent random variables such that Z m i = q if node i of network m belongs to block q with q ∈ {1, . . . , Q m } and

P (Z m i = q) = π m q (1)
where π m q > 0 and Qm q=1 π m q = 1. Given the latent variables Z m , the X m ij 's are assumed to be independent and distributed as

X m ij |Z m i = q, Z m j = r ∼ F(•; α m qr ), ( 2
)
where F is referred to as the emission distribution. F is chosen to be the Bernoulli distribution for binary interactions, and the Poisson distribution for weighted interactions such as counts. Let f be the density of the emission distribution, then:

log f (X m ij ; α m qr ) = X m ij log α m qr + (1 -X m ij ) log 1 -α m qr for Bernoulli emission -α m qr + X m ij log α m qr -log(X ij !)
for Poisson emission .

(3) Equations ( 1), ( 2) and ( 3) define the SBM model and we will now use the following short notation:

X m ∼ F-SBM nm (Q m , π m , α m ). (sep-SBM )
where F encodes the emission distribution, n m is the number of nodes, Q m is the number of blocks in network m, and π m = (π m q ) q=1,...,Qm is the vector of their proportions. The Q m × Q m matrix α m = α m qr q,r=1,...,Qm denotes the connection parameters i.e. the parameters of the emission distribution. Moreover, α m qr ∈ A F where A F = (0, 1) (resp. A F = R * + ) for the Bernoulli (resp. Poisson) emission distribution. In the sep-SBM model, each network m is assumed to follow a SBM with its own parameters (π m , α m ).

First ecological example: three stream food webs We fit the sep-SBM on the 3 stream food webs, respectively referred to as Martins, Cooper and Herlzier. To do so, we use the sbm R-package [START_REF] Chiquet | sbm: Stochastic Blockmodels[END_REF][START_REF] Leger | blockmodels: Latent and Stochastic Block Model Estimation by a 'V-EM' Algorithm. R package version 1[END_REF] on each network, which implements a variational version of the EM algorithm to estimate the parameters and selects the number of blocks Q m using a penalized likelihood criterion ICL. These inference tools (variational EM and ICL) will be recalled hereafter.

We obtain respectively Q 1 = 5 blocks for Martins, Q 2 = 3 blocks for Cooper and Q 3 = 4 blocks for Herlzier. The adjacency matrices of the food webs reordered by block membership are plotted in Figure 1. The two bottom blocks of each food web is composed of basal species (species not feeding on other species). For Cooper, the higher trophic levels are grouped together in the same block: the lack of statistical power does not allow to divide the nodes into more blocks. For Herlzier the higher trophic level is separated into 2 blocks mainly determined on how much they prey on the less preyed basal block. Martins has a separation into 3 blocks, the third one is a medium trophic level, which preys on basal species and is highly preyed on by species of the first block. The first two blocks are made up of higher trophic level species, with the last two blocks being much less connected than the first.

As can be seen in Figure 1, the connectivity structures of these three networks seem to have a lot of similarities. To explore further this aspect, Section 3 is dedicated to the presentation of several colSBM models assuming common structures among the networks of a given collection. A second ecological example: 67 predation networks The three previously presented food webs were extracted from a larger collection which involves more networks collected in similar ecosystems. When willing to analyze a larger collection, with more heterogeneous conditions, one cannot expect to find a structure that will fit all networks well, but it could fit a sub-collection of networks. To illustrate this, we consider a collection of 67 predation networks which are all directed networks with more than 30 species, from the Mangal database [START_REF] Vissault | rmangal: An R package to interact with Mangal database[END_REF]. They are issued from 33 datasets each containing from 1 to 10 networks. The number of species ranges from 31 to 106 (3395 in total) by networks; the networks have density ranging from .01 to .32 (14934 total predation links). Fitting sep-SBM on this collection leads to networks having between 2 to 8 blocks. The number of blocks containing only basal species varies from 0 to 3 betwwen the 67 networks which leads to quite contrasted structures. A method to cluster the networks of this collection will be developed in Section 6.

Joint modeling of a collection of networks

We now present a set of probabilistic models designed to introduce structure consensus into a collection of networks of interest. For ease of notation, we develop the models for directed networks; extensions to the undirected cases are straightforward. Note that the networks (X m ) m=1,...,M are always assumed to be independent random objects. A summary of the various models is provided in Table 1, from the most to the less constrained model

A collection of i.i.d. SBM

The first model we propose is the most constrained one and assumes that the networks are independent realizations of the same Q-block SBM model with identical parameters. The so-called iid-colSBM states that: (iid-colSBM ) where ∀(q, r) ∈ {1, . . . , Q} 2 , α qr ∈ A F , π q ∈ (0, 1] and Q q=1 π q = 1. The model involves (Q-1)+Q 2 parameters, the first term corresponding to the block proportions (π 1 . . . , π Q-1 ) and the second term to the connection parameters.

X m ∼ F-SBM nm (Q, π, α), ∀m = 1, . . . M,
However, assuming that the blocks are represented in the same proportions in each network is a strong assumption that may lead to the model being of little practical use. In food webs, the proportion of species at a given trophic level may differ between networks that nevertheless share the same structure, for example networks coming from various studies may have different species resolution (the number of basal species may differ from a network to another). The following model relaxes this assumption.

A collection of networks with varying block sizes

π-colSBM still assumes that the networks share a common connectivity structure encoded in α, but that the proportions of the blocks are specific to each network. More precisely, for m ∈ {1, . . . , M }, the X m are independent and

X m ∼ F-SBM nm (Q, π m , α) ∀m = 1, . . . M. (π-colSBM )
where ∀(q, r) ∈ {1, . . . , Q} 2 , α qr ∈ A F and Q q=1 π m q = 1, ∀m ∈ {1, . . . , M }. In order to make the model more flexible and suitable to ecological networks, we allow some block proportions π m q to be null in certain networks (π m q ∈ [0, 1]): if π m q = 0 then block q is not represented in network m. The connectivity structure of each network is then part of a larger connectivity structure common to all networks of the collection.

In order to ensure the identifiability of the model, we assume that each block q is represented in at least one network m, or equivalently, for any block q ∈ {1, . . . , Q}, ∃m ∈ {1, . . . , M } such that π m q > 0. Let S be the M × Q support matrix such that ∀(m, q) S mq = 1 π m q >0 . Then, the set of admissible supports is the set of binary matrices with at least one 1 in each row and column:

S Q := S ∈ M M,Q ({0, 1}), M m=1 S mq ≥ 1 ∀q = 1, . . . , Q, Q q=1 S mq ≥ 1 ∀m = 1, . . . , M
For a given number of blocks Q and matrix S, the number of parameters of the π-colSBM model is deduced as follows:

NP(π-colSBM ) = M m=1 Q q=1 S mq -1 + Q q,r=1 1 (S ′ S)qr>0
The first term corresponds to the non-null block proportions in each network. The second quantity accounts for the fact that some blocks may never be represented simultaneously in any network, so the corresponding connection parameters α qr are not useful for defining the model (see the illustration below).

Remark 1. Note that NP(π-colSBM ) ≤ M (Q -1) + Q 2 , this upper bound corresponding to the case where all the blocks are represented in all the networks, but with varying proportions.

Illustration We illustrate the flexibility of π-colSBM model with three examples, all with Q = 3 and M = 2.

1. First consider the situation where the 3 blocks are represented in the two networks but with different block proportions: .20, .50, .30] .

α =   α 11 α 12 α 13 α 21 α 22 α 23 α 31 α 32 α 33   π 1 = [.25, .25, .50] π 2 = [
In that case, S = 1 1 1 1 1 1 and the number of parameters is 2(3-1)+3×3 = 13.

2. Now imagine two networks with nested structures. Blocks 1 and 3 are represented in the two networks while block 2 only exists in network 1. In this illustration, block 2 may refer to a block of parasites which are not always included in food webs [START_REF] Lafferty | Parasites in food webs: the ultimate missing links[END_REF]. In that case, S = 1 1 1 1 0 1 and the number of parameters is (3 -1) + (2 -

1) + 3 × 3 = 12.
3. Finally, let us consider two networks with partially overlapping structures.

The two networks share block 1 (for instance basal species) but the remaining nodes of each network cannot be considered as equivalent in terms of connectivity. One may think of species belonging to trophic chains with different connectivity patterns. .40, 0 , .60] .

α =   α 11 α 12 α 13 α 21 α 22 • α 31 • α 33   π 1 = [.25, .75, 0 ] π 2 = [
In that case, S = 1 1 0 1 0 1 . Moreover, blocks 2 and 3 never interact because their elements do not belong to the same network and so α 23 and α 32 are not required to define the model. As a consequence, the number of parameters is equal to (2 -1) + (2 -1) + 7 = 9.

Model name Block prop.

Connection param. Nb of param. iid-colSBM π m q = π q , π q > 0 α m qr = α qr (Q -1) + Q 2 π-colSBM π m q , π m q ≥ 0 α m qr = α qr ≤ M (Q -1) + Q 2 δ-colSBM π m q = π q , π q > 0 α m qr = δ m α qr (Q -1) + Q 2 + (M -1) δπ-colSBM π m q , π m q ≥ 0 α m qr = δ m α qr ≤ M Q + Q 2 -1 sep-SBM π m q , π m q > 0 α m qr M m=1 (Q m -1) + Q 2 m
Table 1: Summary of the various models defined in Section 3. The last line corresponds to modeling separately each network as presented in Section 2.

A collection of networks with varying density (δ-colSBM )

The iid-colSBM can be relaxed in another direction, assuming that the M networks exhibit similar intra-and inter-blocks connectivity patterns but with different densities. More precisely, let δ m ∈ R be a density parameter for network m. The δ-colSBM is defined as follows:

X m ∼ F-SBM nm (Q, π, δ m α). (δ-colSBM ) with π q > 0, ∀q = 1, . . . , Q, Q q=1 π q = 1. Moreover ∀(m, q, r
), δ m α qr ∈ A F and one of the density parameter equal to one (δ 1 = 1) for identifiability purpose. This model mimics different intensities of connection between networks. In ecology, these differences in densities between networks could be due to different sampling efforts for instance, leading to varying total numbers of observed interactions. δ-colSBM involves NP(δ-colSBM ) = (Q -1) + Q 2 + (M -1) parameters.

Collection of networks with varying block sizes and density (δπ-colSBM )

Finally, we propose to mix the models π-colSBM and δ-colSBM to obtain a more complex one which allows each network to have its own block proportions π m as well as a specific scale density parameter δ m . Then, the (X m ) m∈{1,...,M } are independent and

X m ∼ F-SBM nm (Q, π m , δ m α), (δπ-colSBM ) where ∀(m, q, r), δ m α qr ∈ A F , δ 1 = 1, π m q ≥ 0 and Q q=1 π m q = 1. The number of parameters is given by NP(δπ-colSBM ) = NP(π-colSBM ) + M -1, the last term corresponding to the aditional proper density of each network. Note that NP(δπ-colSBM ) ≤ M (Q -1) + Q 2 + M -1 = M Q + Q 2 -1.

Likelihood and identifiability of the models

In this section, we derive the expression of the likelihood of the most complex model δπ-colSBM and provide conditions to ensure the identifiability of the parameters for each of the four models.

Log-likelihood expression

For a given matrix S, let θ S be:

θ S = (π 1 , . . . π M , δ 1 , . . . , δ M , α) = (π, δ, α),
where π m q = 0 for any q such that S mq = 0. Let Z m iq = 1 Z m i =q be the latent variable such that Z m iq = 1 if node i of network m belongs to block q, Z m iq = 0 otherwise. We define Z m = (Z m iq ) i=1 ...,nm,q=1...,Q . Then the log likelihood is:

ℓ(X; θ S ) = M m=1 log Z m exp {ℓ(X m |Z m ; α, δ) + ℓ(Z m ; π)} dZ m , (4)
where

ℓ(X m |Z m ; α, δ) = nm i,j=1 i̸ =j (q,r)∈Qm Z m iq Z m jr log f X m ij ; δ m α qr , ℓ(Z m ; π) = nm i=1 q∈Qm Z m iq log π m q .
with Q m = {q ∈ {1, . . . , Q}|π m q > 0} and f defined as in Equation ( 3). The log-likelihood functions of the other models can be deduced from this one, setting δ m = 1 for iid-colSBM and π-colSBM and π m = π for iid-colSBM and δ-colSBM with S being a matrix of ones (all blocks are represented in each network).

Identifiability

We aim at giving conditions ensuring the identifiability of the models we propose. We aim at proving that if ℓ(X; θ) = ℓ(X; θ ′ ) for any collection X then θ = θ ′ . The proof relies on the identifiability for the standard SBM demonstrated by [START_REF] Celisse | Consistency of maximum-likelihood and variational estimators in the stochastic block model[END_REF] and is provided in Section 1 of the Supplementary Material [START_REF] Chabert-Liddell | Supplement to "learning common structures in a collection of networks[END_REF]. Note that, like any mixture models, all the models are identifiable up to a label switching of the blocks.

Properties 1.

iid-colSBM The parameters (π, α) are identifiable up to a label switching of the blocks provided that:

(1.1) ∃m * ∈ {1, . . . , M } : n m * ≥ 2Q, (1.2) (α • π) q ̸ = (α • π) r ∀(q, r) ∈ {1, . . . , Q} 2 , q ̸ = r.
δ-colSBM The parameters (π, α, δ 1 , . . . , δ M ) are identifiable up to a label switching of the blocks provided that:

(2.1)

∃m * ∈ {1, . . . , M } : n m * ≥ 2Q and δ m * = 1, (2.2) (α • π) q ̸ = (α • π) r ∀(q, r) ∈ {1, . . . , Q} 2 , q ̸ = r. π-colSBM Assume that ∀m = 1, . . . , M, X m ∼ F-SBM nm (Q, π m , α). Let Q m = |Q m | = |{q = 1 . . . , Q, π m q > 0}|
be the number of non empty blocks in network m. Then the parameters (π 1 , . . . , π M , α) are identifiable up to a label switching of the blocks under the following conditions:

(3.1) ∀m ∈ {1, . . . , M } : n m ≥ 2Q m , (3.2) ∀m ∈ {1, . . . , M }, (α • π m ) q ̸ = (α • π m ) r ∀(q, r) ∈ Q 2 m , q ̸ = r, (3.3) Each diagonal entry of α is unique. δπ-colSBM Assume that ∀m = 1, . . . , M, X m ∼ F-SBM nm (Q, π m , δ m α). Let Q m = |Q m | = |{q = 1 . . . , Q, π m q >
0}| be the number of non empty blocks in network m. Then the parameters (π 1 , . . . , π M , α, δ 1 , . . . , δ M ) are identifiable up to a label switching under the following conditions:

(4.1) ∀m ∈ {1, . . . , M } : n m ≥ 2|Q m |, (4.2) δ 1 = 1, If Q ≥ 2: (4.3) (α • π m ) q ̸ = (α • π m ) r for all (q ̸ = r) ∈ Q 2 m , (4.4) ∀m ∈ {1, . . . , M }, Q m ≥ 2, (4.5) Each diagonal entry of α is unique, If Q ≥ 3:
(4.6) There is no configuration of four indices (q, r, s, t) ∈ {1, . . . , Q} such that α qq /α rr = α ss /α tt with q ̸ = s or r ̸ = t and with q ̸ = r or s ̸ = t,

(4.7) ∀m ≥ 2, |Q m ∩ ∪ l:l<m Q l | ≥ 2.
5 Inference of the models

Variational estimation of the parameters

We now tackle the estimation of the parameters θ S ∈ Θ S for a given support matrix S. For ease of reading, the index S is dropped in this section. The likelihood given in Equation (S-2) is not tractable in practice, even for a small collection of networks as it relies on summing over M m=1 |Q m | nm terms. A well-proven approach to handle this problem for the inference of the SBM is to rely on a variational version of the EM (VEM) algorithm. The approach is similar for both Bernoulli and Poisson models. This is done by maximizing a lower (variational) bound of the log-likelihood of the observed data by approximating p(Z|X; θ) with a distribution on Z named R issued from a family of factorizable distribution [START_REF] Daudin | A mixture model for random graphs[END_REF]:

J (R, θ) := E R [ℓ(X, Z; θ)] + H(R(Z)) ≤ ℓ(X; θ),
where H denotes the entropy of a distribution. The variational distribution R can be fully described by the probabilities τ m iq where

τ m iq = P R (Z m iq = 1). ( 5
)
These quantities approximate the posterior node grouping probabilities.

The VEM algorithm is a two-step iterative procedure which alternates the variational E-step and the M-step. The E-step consists in optimizing J (R, θ) for a current parameter value θ with respect to R constrained to be in the family of factorizable distributions. And the M-step consists in maximizing J (R, θ) with respect to θ for a given variational distribution R. For the colSBM models, networks can be treated independently during the E-step, while the M-step serves as a link between the structures of the networks in the collection. For δ-colSBM and δπ-colSBM when F = Bernoulli, numerical approximations are needed as no explicit expression of δ and α can be derived. Further details of the variational procedure and the expression of the parameters estimators are provided in Section 2 of the Supplementary Material (Chabert-Liddell et al., 2022b).

Model selection

There are two model selection issues. First, under a fixed colSBM , we aim to choose the number of blocks Q and determine the support matrix of the blocks S for π-colSBM and δπ-colSBM . This task is tackled in Subsection 5.2.1 by introducing a penalized likelihood criterion. Second, the comparison of the the colSBM models -each one introducing various degrees of consensus between the networks-with the sep-SBM -which assumes that each network has its own structure-is dealt with in Subsection 5.2.2.

Selecting the number of blocks Q

A classical tool to choose the number of blocks in the SBM context is the Integrated Classified Likelihood (ICL) proposed by [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]; [START_REF] Daudin | A mixture model for random graphs[END_REF]. ICL derives from an asymptotic approximation of the marginal complete likelihood m(X, Z) = θ exp{ℓ(X, Z|θ)}p(θ)dθ where the parameters are integrated out against a prior distribution, resulting in a penalized criterion of the form max θ ℓ(X, Z; θ) -1 2 pen. In the ICL, the latent variables Z are integrated out against an approximation of p(Z|X, θ) obtained via the variational approximation. This leads to the following expression

ICL = max θ E R [ℓ(X, Z; θ)] - 1 2 pen
Using the fact that E R [ℓ(X, Z; θ)] ≈ ℓ(X; θ) -H( R), one understands that, as emphasized in the literature, ICL favors well separated blocks by penalizing for the entropy of the node grouping. However, in this work, our goal is not only to group the nodes into coherent blocks but also to evaluate the similarity of the connectivity patterns between the different networks. As such we would like to authorize models providing grouping of nodes that may be more fuzzy by not penalizing for the entropy. This leads to a BIC-like criterion of the form:

BIC-L = max θ E R [ℓ(X, Z; θ)] + H( R) - 1 2 pen = max θ J ( R, θ) - 1 2 pen
We now supply the expression of the penalty term for the four models we proposed and discuss possible variations of the criterion.

Selection of Q for iid-colSBM and δ-colSBM

For iid-colSBM and δ-colSBM , the derivation of the penalty is a straightforward extension of the classical SBM model, leading to:

BIC-L(X, Q) = max θ J ( R, θ) - 1 2 [pen π (Q) + pen α (Q) + pen δ (Q)] ,
where

pen π (Q) = (Q -1) log M m=1 n m , pen α (Q) = Q 2 log(N M ), pen δ (Q) = 0 for iid-colSBM (M -1) log (N M ) for δ-colSBM .
where

N M = M m=1 n m (n m -1) (6)
is the number of possible interactions. The first term pen π (Q) corresponds to the grouping part where the Q -1 block proportions have to be estimated from the M m=1 n m nodes. The terms pen α (Q) and pen δ (Q) are linked to the connection parameters. Finally, Q is chosen as:

Q = argmax Q∈{1,...,Qmax} BIC-L(X, Q) .

Selection of Q for π-colSBM and δπ-colSBM

Here, in addition to the choice of Q, the collection of support matrices S is considered. In order to penalize the complexity of the model space, we introduce a prior distribution on S defined as follows. Let us introduce Q m = Q q=1 S mq the number of blocks represented in network m. Assuming independent uniform prior distributions on the (Q m )'s and a uniform prior distribution on S for fixed numbers of blocks Q 1 , . . . , Q M represented in each network, we obtain the following prior distribution on S:

log p Q (S) = -M log(Q) - M m=1 log Q Q m
where Q Qm is the number of choices of Q m non-empty blocks among the Q possible blocks in network m. Now, combining the Laplace asymptotic approximation of the marginal complete likelihood (where the parameters have been integrated out) and introducing the prior distribution on S, we obtain the following penalized criterion:

BIC-L(X, Q) = max S max θ S ∈Θ S J ( R, θ S ) - pen π (Q, S) + pen α (Q, S) + pen δ (Q, S) + pen S (Q) 2 ,
where

pen π (Q, S) = M m=1 (Q m -1) log(n m ), pen α (Q, S) = Q q,r=1 1 (S ′ S)qr>0 log (N M ) , pen δ (Q, S) = 0 for π-colSBM (M -1) log (N M ) for δπ-colSBM , pen S (Q) = -2 log p Q (S),
and N M has been defined in Equation ( 6). Finally, Q is chosen as the number of blocks which maximizes the BIC-L criterion:

Q = argmax Q∈{1,...,Qmax} BIC-L(X, Q).
The details about the derivation of this criterion are provided in Section 3 of the Supplementary Material [START_REF] Chabert-Liddell | Supplement to "learning common structures in a collection of networks[END_REF].

Practical model selection

The practical choice of Q and the estimation of its parameters are computationally intensive tasks. Indeed, we should compare all the possible models through the chosen model selection criterion. Furthermore, for each model, the variational EM algorithm should be initialized at a large number of initialization points (due to its sensitivity to the starting point), resulting in an unreasonable computational cost. Instead, we propose to adopt a stepwise strategy, resulting in a faster exploration of the model space, combined with efficient initializations of the variational EM algorithm. The procedure we suggest is given in Algorithm 1 and is implemented in an R-package colSBM available on GitHub: https: //github.com/Chabert-Liddell/colSBM. To initialize a colSBM with Q blocks, we first adjust a sep-SBM with Q blocks, then the Q blocks of the M networks must be associated. This association step can be done in many ways due to label switching within each network which provides us with a lot of possible initializations. Then, the stepwise procedure explores the possible number of blocks by building on the previously fitted models. Note that when fitting the π-colSBM or the δπ-colSBM , the support S has to be determined which is done through an extra-step that consists in thresholding the parameters π m related with the block proportions leading to an exploration over the set S Q .

Algorithm 1: Model selection algorithm Data: X a collection of networks.

begin initialization -Infer sep-SBM on X, with Q ∈ [Q min , Q max ] -Get Ẑm sep-SBM (Q) -Fit
BIC-L is increasing do -Forward loop for Q = Q min + 1, . . . , Q max do -Fit colSBM with Q blocks from initializations obtained by splitting a block in models with Q -1 blocks if π-colSBM or δπ-colSBM then -Fit colSBM with S qm = 1 πm q >t for different value of threshold t -Backward loop for Q = Q max -1, . . . , Q min do
-Fit colSBM with Q blocks from initializations obtained by merging two blocks in models with Q + 1 blocks if π-colSBM or δπ-colSBM then -Fit colSBM with S qm = 1 πm q >t for different value of threshold t -Among all fitted models, keep the b fitted models with the highest BIC-L for each Q return Q = arg max BIC-L(X, Q), with the corresponding θ, Z and S for π-colSBM and δπ-colSBM .

Testing common connectivity structure

We can also use a model selection approach to choose which model from the 4 colSBM s and the sep-SBM is the most adapted to the collection. The most interesting comparison is to decide whether a collection of networks share the same connectivity structure by comparing the model selection criterion obtained for a given colSBM model with the one of sep-SBM . We decide that a collection of networks share the same connectivity structure if:

max Q BIC-L colSBM (X, Q) > M m=1 max Qm BIC-L SBM (X m , Q m ).

Simulation studies

In Section 4 of the Supplementary Material [START_REF] Chabert-Liddell | Supplement to "learning common structures in a collection of networks[END_REF], we perform a large simulation study in order to test the inference and model selection procedures proposed in this section. More specifically, simulating from a π-colSBM model, for various strengths of connectivity structures, we look at our capacity to recover the true connectivity parameter α as well as grouping of the nodes and the true support S. We also assess the quality and the limit of the model selection with the BIC-L criterion. First, we test our ability to select the true number of blocks and to distinguish π-colSBM from iid-colSBM and sep-SBM , and second we vary the block proportions and study its influence on model comparison and the selection of the number of blocks and their support.

We perform another simulation study to understand how, for particular configuration, using a colSBM model on a collection of networks favors the transfer of information between networks and allows to find finer block structures on the networks. An additional one is dedicated to collection of networks simulated from sep-SBM with heterogeneous numbers of nodes. We study how the different colSBM models deal with spurious structures and our ability to detect them through the BIC-L criterion.

6 Partition of a collection of networks according to their connectivity structures

If the networks in a collection do not have the same connectivity structure, we aim to cluster them accordingly. In order to do this, we propose to use the BIC-L criterion in a similar fashion as we did for testing common connectivity structure in Section 5.2.2. We seek the partition of the collection which maximizes a score based on the BIC-L criterion. In this partition, each sub-collection of networks has its own structure given by a colSBM , which represents the best way to model the collection according to the criterion. Clustering a collection of networks consists in finding a partition G = (M g ) g=1,...,G of {1, . . . , M }. Given G, we set the following model on X:

∀g ∈ {1, . . . , G}, ∀m ∈ M g , X m ∼ F-SBM(Q g , π m , δ m α g )
with δ 1 = 1. Moreover, δ m = 1 for all m for iid-colSBM and π-colSBM s and π m = π g for iid-colSBM and δ-colSBM . In other words, the networks belonging to the sub-collection M g share the same mesoscale structure given by a particular colSBM . To any partition G, we associate the following score:

Sc(G) = G g=1 max Q g =1,...,Qmax BIC-L((X m ) m∈Mg , Q g ). ( 7
)
where BIC-L((X m ) m∈Mg , Q g ) is the BIC-L computed on the sub-collection of networks M g . The best partition is chosen as the one which maximizes the score Sc(G) in Equation ( 7).

Computing the BIC-L for all the partitions G requires to consider the 2 M -1 nonempty sub-collections of the networks M, fit the colSBM s on these sub-collections and then combine the associated BIC-L in order to be able to compute the scores given in Equation ( 7). This can be done exhaustively provided that M is not too large but the computational cost becomes prohibitive as M grows.

To circumvent this point, we propose a less computationally intensive forward strategy, starting from G = {1, . . . , M } , and then progressively splitting the collection of networks. In order to explore the space of partitions of {1, . . . , M }, we define a dissimilarity measure between any pair of networks (m, m ′ ) in a subcollection.

This dissimilarity is a squared distance weighted by the block proportions between the connectivity matrices of the two networks. The parameters (block proportions and connectivity matrices) are computed separately on the two networks with the node grouping provided by the inference on the whole sub-collection. We then use 2-medoids clustering to split the sub-collection of networks based on the dissimilarity measures. A split is validated if it increases the score of Equation ( 7). The mathematical definition of the dissimilarity measure and details on the recursive clustering algorithm are given in Appendix A.

Simulation study We illustrate our capacity to perform a partition of a collection of networks based on their structure for all 4 colSBM models in Section 4 of the Supplementary Material [START_REF] Chabert-Liddell | Supplement to "learning common structures in a collection of networks[END_REF].

Application to Food Webs

In this section, we demonstrate the interest of our models on the two collections of ecological networks described in Section 2.

Joint analysis of 3 stream food webs

In Section 2, we fitted sep-SBM and obtained 5 blocks for Martins, 3 blocks for Cooper and 4 blocks for Herlzier. For reminder, a matricial representation of the block reordered food webs was shown in Figure 1. Each food web has 2 blocks of basal species (the 2 bottom blocks).

Finding a common structure between the networks We now fit the four colSBM models in order to find a common structure among the 3 networks. First, notice that our model selection criterion greatly favors common network structure above separated one: BIC-L = -2080 for sep-SBM versus respectively -1964, -1983, -1970 and -1988 for iid-colSBM , π-colSBM , δ-colSBM , δπ-colSBM . The mesoscale structures of the collection under the different models are represented in Figure 2. In this figure, the red shaded matrices represent the estimated connectivity matrices α. In these matrices, the sizes of the blocks are proportional to the block proportion parameters π for iid-colSBM and δ-colSBM and proportional to the averages of block proportion parameters over the networks π m weighted by the number of species within the networks. The cumulative bar plot on the right represents the actual block proportions for each network resulting from the inference of a colSBM denoted as π m where

π m q = nm i=1 τ m iq .
Note that although the iid-colSBM and δ-colSBM assume common parameters for block proportions, the actual block proportions of each network π m fluctuate around the block proportions π .

For each model, the basal species are separated into 2 blocks (bottom blue blocks in Figure 2), similar to the one obtained with the SBMs. For higher trophic levels, the inferred structures slightly differ and are more detailed for the Cooper and Herlzier networks than the ones obtained with SBMs as described hereafter.

• iid-colSBM highlights 5 blocks in total. Block 3 (light green) is a small block of intermediate trophic level species (ones that prey on basal species and are being preyed on by higher trophic levels) with some within block predation. The higher trophic level is divided into 2 more blocks, block 2 (dark green) only preys on the 2 basal blocks, while block 1 (pink) preys on the intermediate block 3 level but only on the most connected basal species block.

• π-colSBM leads also to 5 blocks.There are no empty blocks and the block proportions are roughly corresponding to the ones of iid-colSBM . This renders the flexibility of the π-colSBM of little use compared to the iid-colSBM on this collection.

• With δ-colSBM , the species are grouped into 6 blocks and the networks have different estimated density parameters: δ ≈ (1, .9, 1.2). Block 1 (red), 2 (pink) and 3 (dark green) correspond to the top trophic levels. Block 4 (light green) is an intermediate trophic level group, well connected with both block 1 and the basal species blocks. Block 2 (pink) is huge and only preys on block 6 while block 3 (dark green) is a small group of species that preys on both basal species blocks.

• Finally, δπ-colSBM groups the species into 5 blocks and the networks have different estimated density parameters: δ ≈ (1, .9, 1.4). The connectivity structure is almost similar to the one of π-colSBM but the blocks have different proportions. Block 1 (pink) corresponds to block 1 (red) and 2 (pink) of δ-colSBM and is the merge of two top trophic levels. Again, the block proportions are still quite homogeneous between networks, and the added flexibility of δπ-colSBM compared to δ-colSBM is not used.

Remark 2. On this collection, the entropy of the block memberships is much lower on δπ-colSBM than on δ-colSBM . In this model, to ensure homogeneous block proportions between networks, some nodes tend to get a fuzzy grouping and sit between several blocks (the variational parameters do not concentrate on one block). This phenomenon is taken into account by our model selection criterion which tends to favor models with higher entropy than models with well separated blocks.

Prediction of missing links and dyads Since we have been able to find some common structures between the 3 networks, we now examine if these structures could be used to help recover some information on networks with incomplete information. We proceed as follows: we choose a network m and remove a proportion The ordering is done by trophic level from bottom to top and right to left. For (δ-δπ)colSBM s we give δ below the barplot.

K ∈ [.1, .8] of
• the existing links uniformly at random for the missing link experiment

• or of the existing dyads (both 0 and 1) by encoding them as NA for the missing dyads experiment.

Then, for the missing link experiment, we try to recover where the missing links are among all non existing ones. For the missing dyad experiment, we predict the probability of existence of missing dyads (NA entries). Under the colSBM , the probability of a link between species i and j for network m is predicted by:

p m ij = q,r∈ Qm τ m iq τ m jr δ m α qr ,
where τ m iq and τ m jr are defined as in Section S-2. We resort to the area under the ROC curve to evaluate the capacity of the different models to recover this information. For each value of K, each experiment is repeated 30 times and the results are shown in Figure 3.

First, let us notice that these stream food webs networks have a structure that is well explained by an SBM. When there is little information missing (K < .3) the ROC AUC is over 0.9. Besides, with 70% of missing links or dyads, we still predict better than a random guess (ROC AUC ≈ .75). As there is a common structure between the 3 networks, there is a lot of information to be taken from the ones with no missing information. Starting from K ≥ .2, the colSBM s outperform the sep-SBM on both experiments. Even for K = .8, the prediction is still high. About the difference between the colSBM s, for the missing links experiment, as we remove links from one of the network, its density decreases and the models with a density parameters, δ-colSBM and δπ-colSBM have a built-in mechanism that compensates this fact. As a consequence, these models yield to better predictions than iid-colSBM or π-colSBM for large values of K.

Another noteworthy comment on both the NA and missing links experiments is that as K grows, the amount of information on the modified network gets lower. Hence, π-colSBM and δπ-colSBM lack the statistical power to separate blocks and will empty blocks on this network. This affects our capability to predict the trophic links. On the other hand, iid-colSBM and δ-colSBM will force some separation of the species into blocks, and as the information from the other networks is relevant it still has good prediction performance for large K.

Partition of a collection of 67 predation networks

Now, we consider the collection of networks issued from the Mangal database (Vissault et al., 2020) introduced at page 7 of Section 2. This dataset is too heterogeneous to find a common structure that will fit well on all the networks. Therefore we propose to use a π-colSBM to look for a partition of the networks into groups sharing common connectivity structures. We focus our investigations on this model since we aim to cluster together networks which share some blocks with similar features but we do not expect all the networks in a sub-collection to share exactly the same blocks.

Fitting a π-colSBM on the whole collection provides a 13 blocks connectivity structure with more than half of the blocks being empty (457/871). It leads to a much higher BIC-L than the one given by sep-SBM (-31303 vs. -33311), still the partition we provide below greatly improves this criterion (-30703). In this partition, the networks are clustered into 5 sub-collections. The obtained partition and the connectivity structure of each sub-collection are shown in Figure 4 as well as a contingency table of the obtained sub-collections crossed with the different datasets of the Mangal database. The number of blocks of each sub-collection varies between 8 and 12. The sub-collection denoted A in the following has 12 blocks and contains networks that are denser (networks density ranging from 0.17 to 0.32) than the networks of the other sub-collections (networks density ranging from 0.01 to 0.17). Each sub-collection contains between 1 and 3 blocks which can be considered as blocks containing mainly basal species. Indeed, those blocks have a very low in-going interaction probability with all the other blocks (inferior to 0.02). All the structures exhibit low within blocks connectivity for most blocks ( α g qq is small for most q), meaning that predation links between species of the same block are unlikely, with the exception of the sub-collection A. Also, the networks in all the sub-collections contain mostly trophic chains and only a few cycles, since each of the connectivity parameter matrices α g can almost be reordered as triangular matrices. Our detailed comments on each sub-collection follow.

A This sub-collection consists of 7 networks and 12 blocks are required to describe this sub-collection. 5 networks are issued from the same dataset (id: 80). These 5 networks populate the 12 blocks, while the other 2 networks only populate parts of them. The average density is about 0.18. From the ecological point of view, the blocks can be divided into 3 heterogeneous sets: block 1 to 3 represent the higher trophic levels, block 4 to 8 the intermediate ones and block 9 to 12 the lower ones.

B A sub-collection of 26 networks with heterogeneous size and density issued from various datasets. Most networks populate only parts of the 8 blocks.

The structure is mainly guided by blocks 2, 3, 5 and 6, from higher to lower trophic levels. Block 4 is represented in only 5 networks where it is either an intermediate or a bottom trophic level. It introduces some symmetry in the connectivity matrix rendering it difficult to order the blocks by trophic order. Species from top trophic levels prey on basal species.

C A small sub-collection of 6 networks with density ranging from .06 to .11. All networks are represented in 5 or 6 of the 7 blocks, including the first three blocks. The sub-collection consists of 3 of the 5 networks of dataset 48, the separation being based on the collecting sites. (blocks 6, 7 and 8), while species from block 5 just prey on species from the last two blocks.

E The last sub-collection consists of 6 networks from various datasets. In the 7 blocks structure, the species of block 1 (represented on 4 of the 6 networks) prey on species from all other blocks with the exception of block 7. The basal species are separated between blocks 6 and 7 depending on whether or not they are preyed on by species from the first two blocks.

Discussion

In this paper, we proposed a new method to find a common structure and compare different structures of a collection of networks which we do not assume to share common nodes. This method is very general and could be applied to networks sharing common nodes as well, such as temporal or multiplex networks. Starting from our most basic model, -namely the iid-colSBM -we refined it by proposing models allowing for different mixture distributions and even empty blocks (π-colSBM ), models allowing to find common structure for networks of different density (δ-colSBM ) or even models allowing both (δπ-colSBM ). These are only a few of the possible models within the colSBM framework. Other systems of constraints on the SBM parameters can be proposed. For instance, when driven by the analysis of common community structures, one may impose the diagonal of the connectivity parameter matrix (which corresponds to the within block connectivities) to be common for all the networks but release the restriction on the off-diagonal parameters (between blocks connectivities). The model selection criterion we derived can be used to select the number of blocks but also to choose which colSBM fits better the data. We also presented a strategy providing a partition of a collection of networks into sub-collections of networks sharing common connectivity patterns. The idea behind these models is very general and could be extended to other types of networks. In ecology, bipartite and multipartite networks are common and the model extension [START_REF] Govaert | Clustering with block mixture models[END_REF][START_REF] Bar-Hen | Block models for generalized multipartite networks: Applications in ecology and ethnobiology[END_REF] is straightforward (although some additional modeling choices arise when considering π-colSBM , δ-colSBM or δπ-colSBM ), the main difficulty would then lie in the algorithmic part. Indeed, the dimension of the model space is larger when considering blockmodels for multipartite networks: N * K instead of N * as a number of blocks has to be determined for each of the K different types of units. Thus, our model selection algorithm will not scale well. An adaptation of [START_REF] Côme | Hierarchical clustering with discrete latent variable models and the integrated classification likelihood[END_REF] to the colSBM with the derivation of an exact ICL criterion and a genetic algorithm to search in the model space could be a direction to solve this issue. Additionally, incorporating the type of ecological interaction as a network covariate [START_REF] Mariadassou | Uncovering latent structure in valued graphs: a variational approach[END_REF] would help us understand its impact on the structure of the networks and to the robustness of the ecosystems they depict (Chabert-Liddell et al., 2022a). The main idea of this article could also be extended to the Degree Corrected SBM [START_REF] Karrer | Stochastic blockmodels and community structure in networks[END_REF] which is quite used in practice.

Finally, we notice during our simulations and applications that the colSBM s allow to find a larger number of blocks compared to the sep-SBM and so lead to a finer resolution of the mesoscale structure of the networks. This resolution limit problem was one of the motivations of [START_REF] Peixoto | Hierarchical block structures and high-resolution model selection in large networks[END_REF] and we believe that this direction should be explored further for collections of networks. Also, sometimes practitioners are not interested in the model selection procedures but seek the structure and nodes grouping for a set number of blocks. In the case where networks strongly vary in sizes or reliability of information, it might be interesting to use a weighted likelihood as Škulj and Žiberna (2022) recently proposed for multipartite SBM. Research on how how to choose and cheaply evaluate the weights in the SBM framework remains to be done.

A Details of the network partitioning algorithm

Definition of a dissimilarity measure between networks of a collection (X m ) m∈M It relies on the following steps.

1. Infer colSBM on (X m ) m∈M to get coherent groupings of the nodes encoded in τ m = ( τ m iq ) i=1,...,nm,q=1,..., Q , for any m ∈ M, where τ m iq is defined as in Section S-2. This step supplies a grouping of the nodes into blocks for each network.Note that the inference also supplies the ( δ m ) m∈M , these quantities being set to 1 if we work with the π-colSBM and iid-colSBM .

2. For each network m, compute:

α m qr = nm i,j=1 i̸ =j τ m iq τ m jr X m ij nm i,j=1 i̸ =j τ m iq τ m jr , π m q = nm i=1 τ m iq n m (S-8)
with the convention that π m q = 0 if q / ∈ Q m and α m qr = 0 if {q, r} ̸ ⊂ Q m . These quantities are the separated estimates of the parameters encoding the mesoscale structure for each network, computed from nodes grouping obtained by considering all the networks jointly. They correspond to the parameters estimates of the SBM when the block memberships are known. pi m q corresponds to the expected proportion of nodes of network m grouped in block q, while α m qr is the expected connectivity parameter for network m between group q and r.

3. Then, for any pair of networks (m, m ′ ) ∈ M 2 compute the dissimilarity:

D M (m, m ′ ) = Q (q,r)=1 max π m q , π m ′ q max π m r , π m ′ r α m qr δ m - α m ′ qr δ m ′ 2 . (S-9)
The dissimilarity measure quantifies to what extent the connectivity parameters inferred separately on each network of the pair are different even though the nodes grouping were inferred jointly. This is weighted by the size of the blocks and corrected in the case of δ-colSBM and δπ-colSBM by the density parameter. We take the maximum of the proportion of each block over the two networks in order to further increase the dissimilarity of networks not populating the same blocks. Two networks with low dissimilarity measure between them are expected to have more similar mesoscale structure, and so are more expected to be part of the same sub-collection than two network with higher dissimilarity measure between them.

An algorithm to cluster the collection of networks Now, we use this dissimilarity to guide the search for the best partition of the collection of networks by using Algorithm 2 which consists in a recursive partitioning of the collection. It relies on a 2-medoids algorithm on the dissimilarity matrix defined in Equation (S-9) to obtain a partition of the collection into two sub-collections of networks. This step is repeated recursively until the score, based on the BIC-L, defined in Equation ( 7) calculated from the new partition stops increasing.

Algorithm 2: Clustering a collection of networks into two sub-collections Call: Clust2Coll(X = (X m ) m∈M ) Data: X = (X m ) m∈M a collection of networks and G = {M} the trivial partition in a unique sub-collection begin -Fit colSBM on X -Compute the score Sc 0 = Sc(G) -Compute the dissimilarity for all the networks in the collection

D M (m, m ′ ) m,m ′ ∈M -Apply a 2-medoids algorithm to obtain G 1 and G 2 giving a partition of M. -Compute Sc * = Sc(G * ) where G * = {G 1 , G 2 }. if Sc 0 > Sc * then return G else return Clust2Coll (X m ) m∈G 1 ,Clust2Coll (X m ) m∈G 2

Supplementary Material

S-1 Proof of the identifiability of colSBMs Properties 1.

iid-colSBM The parameters (π, α) are identifiable up to a label switching of the blocks provided that:

(1.1) ∃m * ∈ {1, . . . , M } :

n m * ≥ 2Q, (1.2) (α • π) q ̸ = (α • π) r ∀(q, r) ∈ {1, . . . , Q} 2 , q ̸ = r.
δ-colSBM The parameters (π, α, δ 1 , . . . , δ M ) are identifiable up to a label switching of the blocks provided that:

(2.1)

∃m * ∈ {1, . . . , M } : n m * ≥ 2Q and δ m * = 1, (2.2) (α • π) q ̸ = (α • π) r ∀(q, r) ∈ {1, . . . , Q} 2 , q ̸ = r. π-colSBM Assume that ∀m = 1, . . . , M, X m ∼ F-SBM nm (Q, π m , α). Let Q m = |Q m | = |{q = 1 . . . , Q, π m q > 0}|
be the number of non empty blocks in network m. Then the parameters (π 1 , . . . , π M , α) are identifiable up to a label switching of the blocks under the following conditions:

(3.1) ∀m ∈ {1, . . . , M } : n m ≥ 2Q m , (3.2) ∀m ∈ {1, . . . , M }, (α • π m ) q ̸ = (α • π m ) r ∀(q, r) ∈ Q 2 m , q ̸ = r, (3.3) Each diagonal entry of α is unique. δπ-colSBM Assume that ∀m = 1, . . . , M, X m ∼ F-SBM nm (Q, π m , δ m α). Let Q m = |Q m | = |{q = 1 . . . , Q, π m q > 0}|
be the number of non empty blocks in network m. Then the parameters (π 1 , . . . , π M , α, δ 1 , . . . , δ M ) are identifiable up to a label switching of the blocks under the following conditions:

(4.1) ∀m ∈ {1, . . . , M } : n m ≥ 2|Q m |, (4.2) δ 1 = 1, If Q ≥ 2: (4.3) (α • π m ) q ̸ = (α • π m ) r for all (q ̸ = r) ∈ Q 2 m , (4.4) ∀m ∈ {1, . . . , M }, Q m ≥ 2, (4.5) Each diagonal entry of α is unique, If Q ≥ 3:
(4.6) There is no configuration of four indices (q, r, s, t) ∈ {1, . . . , Q} such that α qq /α rr = α ss /α tt with q ̸ = s or r ̸ = t and with q ̸ = r or s ̸ = t,

(4.7) ∀m ≥ 2, |Q m ∩ ∪ l:l<m Q l | ≥ 2.
Proof. [START_REF] Celisse | Consistency of maximum-likelihood and variational estimators in the stochastic block model[END_REF] proved that the parameters (π m , α m ) of the F-SBM nm (Q m , π m , α m ) are identifiable up to a label switching of the blocks from the observation of a single network X m when F is the Bernoulli distribution and provided that

1. n m > 2Q m , 2. (α m • π m ) q ̸ = (α m • π m ) r for all (q ̸ = r) ∈ {1, . . . , Q m } 2 .
Although they only consider the case where the emission distribution is the Bernoulli distribution, the extension to the Poisson distribution is straightforward. We prove the identifiability of our colSBM models by using their result, F being either the Bernoulli or the Poisson distribution. The proofs for iid-colSBM and δ-colSBM are straightforward while the proofs for π-colSBM and δπ-colSBM are more complicated due to the possible existence of empty clusters in some networks X m .

iid-colSBM Under this model, for all m = 1, . . . M , X m ∼ F-SBM nm (Q, π, α). As a consequence, following [START_REF] Celisse | Consistency of maximum-likelihood and variational estimators in the stochastic block model[END_REF], the identifiability of α and π is derived from the distribution of X m * under assumptions (1.1) and (1.2).

δ-colSBM Under this model, for all m = 1 . . . , M , X m ∼ F-SBM nm (Q, π, δ m α).

Under assumptions (2.1) and (2.2), we obtain the identifiability of α and π from network X m * [START_REF] Celisse | Consistency of maximum-likelihood and variational estimators in the stochastic block model[END_REF]. Now, for any m ̸ = m * , by definition of F-SBM nm (Q, π, δ m α), we have:

E[X m ij ] = δ m π ′ απ .
This proves that δ m is identifiable.

π-colSBM Note that under π-colSBM , we have, for all m,

X m ∼ F-SBM nm (Q m , π m , α m )
where π m is a vector of non-zero proportions of length Q m and α m is the restriction of α to Q m . Under assumptions (3.1) and (3.2) and applying [START_REF] Celisse | Consistency of maximum-likelihood and variational estimators in the stochastic block model[END_REF] we obtain the identifiability of the parameters π m and α m from each network X m separately. However, the identifiability of each π m and α m is established up to a label switching of the blocks in each network. We now have to find a coherent reordering between the networks which takes into account that some blocks are not represented in all the networks. Let us build the complete matrix α using the α m . We fill the diagonal of α which is composed of (diag( α m )) m=1,...,M , taking the unique values and sorting them by increasing order, i.e. α 11 < α 22 < • • • < α QQ . This task is possible because of assumption (3.3). Now, we get back to π m and reorganize them to match with α. For any m, we define ϕ m : {1, . . . , Q m } → {1, . . . , Q} such that α ϕm(q),ϕm(q) = α m qq . With the (ϕ m ) we are able to fill the rest of α as: α ϕm(q)ϕm(r) = α m qr for all (q, r) ∈ {1, . . . , Q m } 2 . Finally, we define π m a vector of size Q such that:

π m q = 0 ∀q ∈ {1, . . . , Q}\ϕ m ({1, . . . , Q m }) π m ϕ -1 m (q) ∀q ∈ ϕ m ({1, . . . , Q m })
. δπ-colSBM We now consider the model where

X m ∼ F-SBM nm (Q, π m , δ m α).
(S-1)

Like in model π-colSBM , Equation (S-1) implies that marginally,

X m ∼ F-SBM nm (Q m , π m , α m
) where α m = δ m (α qr ) q,r∈Qm .

Applying [START_REF] Celisse | Consistency of maximum-likelihood and variational estimators in the stochastic block model[END_REF] on the distribution of X m , we obtain the identifiability of the parameters π m and α m (assumptions (4.1) and (4.3)). We now have to match the structures of the networks and to take into account the empty blocks. We separate the cases where Q = 2 from the ones where Q > 3.

The proof relies on a sequential identification the parameters (with respect to m).

• For Q = 2, we do not allow empty blocks (assumption (4.4), Q m ≥ 2) so α m = δ m α and π m = π m . Using the fact that δ 1 = 1 (assumption (4.2)), we identify π 1 and α. Since we know that the diagonal entries of α are unique (assumption (4.5)), α can be chosen such that α 11 > α 22 . This provides the ordering of the blocks in each network. Then, we identify the π m in a unique manner and not up to label switching.

• For Q ≥ 3, for each m ∈ {1, . . . , M }, by assumptions (4.1) and (4.3) and using the marginal distributions, we are able to identify α m and π m . Using the fact that δ 1 = 1 (assumption (4.2)) and the fact that the entries of the diagonal of α are unique, we can do as in π-colSBM and identify π 1 and (α qr ) (q,r)∈Q 1 .

Then for m = 2, up to a relabelling of the blocks in α 2 , we can define δ 2 = α 2 11 /α 11 = α 2 22 /α 22 since there are at least two blocks in network m = 2 that correspond to two blocks already identified in network m = 1 by assumption (4.7). We then need to prove that this parameter δ 2 is uniquely defined. Let us assume that there exists a parameter δ ′ 2 such that δ ′ 2 ̸ = δ 2 and δ ′ 2 = α 2 ii /α kk = α 2 jj /α ll with i ̸ = j. By definition of the δπ-colSBM α uu := α 2 ii /δ 2 which is either an already identified parameter from network m = 1 or corresponds to a new block represented in network m = 2 but not in network m = 1. Note that, since δ ′ 2 ̸ = δ 2 , then u ̸ = k. For the same reason, ∃v ̸ = l such that α vv = α 2 jj /δ 2 . Since i ̸ = j, the parameters α 2

ii and α 2 jj are not equal which implies that α kk ̸ = α ll and α uu ̸ = α vv and finally that k ̸ = l and u ̸ = v. By computing:

α uu α vv = α 2 ii /δ 2 α 2 jj /δ 2 • δ ′ 2 δ ′ 2 = α kk α ll
we obtain a contradiction with assumption (4.6). Therefore, δ 2 = δ ′ 2 . We can then identify the blocks in network m = 2 by matching α 2 qq /δ 2 with the α qq already identified. The α 2 qq /δ 2 that do not match with the previously identified parameters complete the matrix α. The process is iterated with networks m = 3, . . . , M . Once the matrix α and the parameters in δ are identified, injections from {1, . . . , Q m } → {1, . . . , Q}, corresponding to the matching of the blocks, provide the π m .

S-2 Variational estimation of the parameters

We provide here some details for the estimation of the parameters θ S ∈ Θ S for a given support matrix S. For ease of reading, the index S is dropped in this section. The likelihood

ℓ(X; θ) = M m=1 log Z m exp {ℓ(X m |Z m ; α, δ) + ℓ(Z m ; π)} dZ m (S-2)
is not tractable in practice, even for a small collection of networks as it relies on summing over M m=1 |Q m | nm terms. A well-proven approach to handle this problem for the inference of the SBM is to rely on a variational version of the EM algorithm. This is done by maximizing a lower (variational) bound of the log-likelihood of the observed data [START_REF] Daudin | A mixture model for random graphs[END_REF]. More precisely,

ℓ(X; θ) = M m=1 ℓ(X m ; θ) ≥ M m=1 ℓ(X m ; θ) -D KL (R m (Z m )∥p(Z m |X m ; θ))
where D KL is the Kullback-Leibler divergence and R m stands for any distribution on Z m . The last equation can be reformulated as:

ℓ(X; θ) ≥ M m=1 E Rm [ℓ(X m , Z m ; θ)] + H(R m ) =: J (R, θ).
(S-3)

where R = ⊗ M m=1 R m and H denotes the entropy of a distribution. Now, if, for all m ∈ {1, . . . , M }, R m is chosen in the set of fully factorizable distributions and if one sets τ m iq = P Rm (Z m iq = 1) then H(R m ) is equal to:

H(R m ) = - nm i=1 q∈Qm
τ m iq log τ m iq .

(S-4)

Besides, the complete likelihood of network m for the δπ-colSBM marginalized over R m is given by:

E Rm [ℓ(X m , Z m ; θ)] = nm i,j=1 i̸ =j (q,r)∈Qm τ m iq τ m jr log f (X m ij ; δ m α qr ) + nm i=1 q∈Qm
τ m iq log π m q .

(S-5) Finally, the variational lower bound J (R, θ) := J (τ , θ) is obtained by plugging Equations (S-4) and (S-5) into the right member of Equation (S-3). Note that the lower bound J (τ , θ) is equal to the log-likelihood if R m (Z m ) = p(Z m |X m ; θ) for all m ∈ {1, . . . , M }.

The variational EM (VEM) algorithm consists in optimizing the lower bound J (τ , θ) with respect to (τ , θ), by iterating two optimization steps with respect to τ and θ respectively, also referred to as VE-step and M-step. The details of each step are specific to the model at stake and are detailed hereafter. VE-step At iteration (t) of the VEM algorithm, the VE-step consists in maximizing the lower bound with respect to τ :

τ (t+1) = arg max τ J (τ , θ (t)
).

Note that by doing so, one minimizes the Kullback-Leibler divergences between R m (Z m ) and p(Z m |X m ), and so approximates the true conditional distribution p(Z m |X m ) in the space of fully factorizable probability distributions. The τ m 's can be optimized separately by iterating the following fixed point systems for all m ∈ {1, . . . , M }:

τ m(t+1) iq ∝ π m(t) q nm j=1 j̸ =i r∈Qm f (X m ij ; δ (t) m α (t) qr ) τ m(t+1) jr ∀i = 1, . . . , n m , q ∈ Q m . (S-6) M-
Step At iteration (t) of the VEM algorithm, the M-step maximizes the variational bound with respect to the model parameters θ:

θ (t+1) = arg max θ J ( τ (t+1) , θ).
The update depends on the chosen model and the estimations are derived by canceling the gradient of the lower bound. For the sake of simplicity, the iteration index (t) is dropped in the following formulae. The obtained formulae involve the following quantities:

e m qr = nm i,j=1 i̸ =j τ m iq τ m jr X m ij , n m qr = nm i,j=1 i̸ =j τ m iq τ m jr , n m q = nm i=1 τ m iq .
On the one hand, the (π

(m)
q ) q∈Qm are estimated as

π m q =
n m q n m for π-colSBM and δπ-colSBM , which is the expected proportion of the nodes in each allowed block for network m.

On the other hand,

π q = M m=1 n m q M m=1 n m
for iid-colSBM and δ-colSBM , taking into account all the networks at the same time. The connection parameters α qr of iid-colSBM and π-colSBM are estimated as the ratio of the number of interactions between blocks q and r among all networks over the number of possible interactions:

α qr = M m=1 e m qr M
m=1 n m qr for iid-colSBM and π-colSBM .

For the δ-colSBM and δπ-colSBM , there is no closed form for δ and α for a given value of τ . If F = Poisson, then δ and α can be iteratively updated using the following formulae: Remark S-1. In practice when F = Bernoulli, to update the parameters δ and α we use by default the method of moving asymptotes [START_REF] Svanberg | A class of globally convergent optimization methods based on conservative convex separable approximations[END_REF]. A much faster method consists in making a non constrained optimization by using the same estimates as the ones for F = Poisson. The estimates α qr δ m are clipped to (0, 1) for all q, r, m. Users of the colSBM package can choose which of these two methods to use as well as any of the ones implemented in the NLOPT package [START_REF] Johnson | The nlopt nonlinear-optimization package[END_REF].

Remark S-2. In the VE-Step, each network can be treated independently, so the computation can be parallelized with ease. Also, it can be more efficient to update only a subset of networks at each step to avoid being stuck in local maxima. So we use a slightly modified VEM algorithm where we just compute the VE-step on one network at a time (the order of which is taken uniformly at random) before updating the corresponding parameters in the M-Step.

S-3 Details of the model selection procedure when allowing for empty blocks

For π-colSBM and δπ-colSBM , the model is described by its support S. We can compute the likelihood for a given support. We recall that θ S = {α S , δ, π S } are the parameters restricted to their support. Then for the model represented by S, the complete likelihood is given by:

p(X, Z|S) = θ S p(X, Z|θ S , S)p(θ S )d(θ S ) = (α S ,δ) π S p(X|Z, α S , δ, S)p(Z|π S , S)p(α S , δ)p(π S )d(α S , δ)d(π S ) = (α S ,δ) p(X|Z, α S , δ, S)p(α S , δ)d(α S , δ) B1 π S p(Z|π S , S)p(π S )d(π S ) B2 , (S-7)
where the prior on the emission parameters and on the mixture parameters are assumed to be independent. The restriction of the parameter space to the one associated with the support S is needed. Otherwise, some parameters would not be defined or would lie on the boundary of the parameters space, and the following asymptotic derivation would not be properly defined. We use a BIC approximation on B1 where we rewrite:

p(X|Z, S) = (α S ,δ) M m=1 p(X m |Z m , α S , δ, S) p(α S , δ)d(α S , δ) = max (α S ,δ) exp M m=1 ℓ(X m |Z m ; α S , δ, S) - 1 2 ν(α S , δ) log m n m (n m -1) + O(1) , where ν(α S , δ) = ν(α S ) = Q q,r=1 1 (S ′ S)qr>0 , for π-colSBM ν(α S ) + M -1 for δπ-colSBM.
For B2, we use a Q m -dimensional Dirichlet prior for each π m S :

p(Z|S) = m∈M π m S p(Z m |π m S )p(π m S )d(π m S ) = max π S exp M m=1 ℓ(Z m ; π m S ) - Q m -1 2 log(n m ) + O(1) .
Then, we input B1 and B2 into Equation (S-7) to obtain:

log p(X, Z|S) ≈ max θ S ℓ(X, Z; θ|S)- 1 2 pen π (Q, S)+pen α (Q, S)+pen δ (Q, S) , (S-8)
with the penalty terms as given in the main text.

As Z is unknown we replace each Z m iq by the variational parameters τ m iq which maximizes the variational bound for a given support S. Then, we add the entropy of the variational distribution H( R) to Equation (S-8). This leads to the variational bound of Equation (S-3), as

max θ S J (τ , θ S ) = max θ S ℓ(X, E R [Z]; θ S ) + H( R),
which we recall is a surrogate of the log-likelihood of the observed data. We define

BIC-L(X, Q, S) = max θ S J (τ , θ|S) - 1 2 pen π (Q, S) + pen α (Q, S) + pen δ (Q, S)
which is a penalized likelihood criterion when the support S is known.

Finally to obtain the criterion BIC-L(X, Q) for π-colSBM and δπ-colSBM , we need to penalize for the size of the space of possible models that depends on the support S. For a given Q corresponding to the number of different blocks in the collection of networks X, we set the prior on S decomposed as the product of uniform priors on the numbers of blocks (between 1 and Q) actually represented in each network and uniform priors for the choice of these Q m blocks among the Q possible blocks (Q m is the number of blocks that are represented in network m):

p Q (S) = p Q (Q 1 , . . . , Q M ) • p Q (S|Q 1 , . . . , Q M ) = 1 Q M • M m=1 1 Q Q m .
The prior is given on the space S Q of admissible support. Using a BIC approximation and under a concentration assumption on the correct support, we derive

log p(X, Z|Q) = log S p(X, Z|S, Q)p Q (S)dS ≈ log S exp BIC-L(X, Q, S)p Q (S) dS ≈ max S∈S Q BIC-L(X, Q, S) -log p Q (S) .
Thus, by denoting pen S (Q) = -2 log p Q (S) in the equation above, we obtain:

BIC-L(X, Q) = max S∈S Q max θ S J ( τ , θ S ) - 1 2 [pen π (Q, S) + pen α (Q, S) + pen δ (Q, S) + pen S (Q)] ,
where

pen π (Q, S) = M m=1 (Q m -1) log(n m ), pen α (Q, S) = Q q,r=1 1 (S ′ S)qr>0 log (N M ) , pen δ (Q, S) = 0 for π-colSBM (M -1) log (N M ) for δπ-colSBM , pen S (Q) = -2 log p Q (S).

S-4 Simulation studies

In this section, we perform a large simulation study. The first study aims at testing the ability of the inference method to recover the number of blocks and the parameters for the π-colSBM model. The second study highlights the performances in terms of clustering of networks based on their mesoscale structure.

S-4.1 Efficiency of the inference procedure

Simulation paradigm Let us simulate data under the π-colSBM model with M = 2, n m = 120 and Q = 4. α and π are chosen as:

α = .25 +     3ϵ α 2ϵ α ϵ α -ϵ α 2ϵ α 2ϵ α -ϵ α ϵ α ϵ α -ϵ α ϵ α 2ϵ α -ϵ α ϵ α 2ϵ α 0     , π 1 = σ 1 (.2, .4, .4, 0), π 2 = σ 2 (0, 1 3 , 1 3 , 1 3 ).
(S-9) with ϵ α taking eight equally spaced values ranging from 0 to 0.24. For each value of ϵ α , 30 datasets (X 1 , X 2 ) are simulated, resulting in 8 × 30 = 240 datasets. More precisely, for each dataset, we pick uniformly at random two permutations of {1, . . . , 4} (σ 1 , σ 2 ) with the constraint that σ 1 (4) ̸ = σ 2 (1). This ensures that each of the two networks have a non-empty block that is empty in the other one. Then the networks are simulated with Bern-SBM 120 (4, α, π m ) with the previous parameters.

Each network has 2 blocks in common and their connectivity structures encompass a mix of core-periphery, assortative community and disassortative community structures, depending on which 3 of the 4 blocks are selected for each network. ϵ α represents the strength of these structures, the larger, the easier it is to tell apart one block from another.

Inference On each simulated dataset, we fit the iid-colSBM , π-colSBM and sep-SBM models. The inference is performed with the VEM algorithm and the BIC-L criterions presented in the main manuscript.

Quality indicators

The assess the quality of the inference, we compute the following set of indicators for each simulated dataset.

• First, for each dataset, we put in competition π-colSBM with sep-SBM and iid-colSBM respectively. To do so, for each dataset, we compute the BIC-L of each model π-colSBM is preferred to sep-SBM (resp. iid-colSBM ) if its BIC-L is greater.

• Secondly, when considering the π-colSBM , we compare Q to its true value (Q = 4).

• For π-colSBM and Q fixed to its true value (Q = 4), we evaluate the quality of recovery of the support matrix S by calculating:

Rec( S, S) = max σ∈S 4

1 {∀q,mSmq= S mσ(q) } (S-10) the greater the better.

• In order to evaluate the ability to recover the true connectivity parameter in the π-colSBM model, we compare α to its true value for the true number of blocks Q = 4 through: RMSE( α, α) = min ( α σ(q)σ(r) -α qr ) 2 , (S-11) the σ being there to correct the possible label switching of the blocks.

• Finally, we judge the quality of our grouping of the nodes into blocks with the Adjusted Rand Index [START_REF] Hubert | Comparing partitions[END_REF], ARI = 0 for a random grouping and 1 for a perfect recovery). For each network, for the π-colSBM , using Q, we compare the block memberships to the real ones by taking the average over the two networks

ARI = 1 2 ARI( Z 1 , Z 1 ) + ARI( Z 2 , Z 2 )
and by computing it on the whole set of nodes

ARI 1,2 = ARI ( Z 1 , Z 2
), (Z 1 , Z 2 ) .

All these quality indicators are averaged among the 30 simulated datasets. The results are provided in Table S-1. Each line corresponds to the 30 datasets simulated with a given value of ϵ α . The first columns concatenate the results of the model comparison task. The following set of columns is about the selection of Q and the estimation of S. The last columns supply the RMSE on α and the ARI.

where σ m is a permutation of {1, 2, 3} proper to network m and σ(π) = (π σ(i) ) i=1,...,3 . The networks are divided into 3 sub-collections of 3 networks with connectivity parameters as follows: .4]. α as represents a classical assortative community structure, while α cp is a layered core-periphery structure with block 2 acting as a semi-core. Finally, dis is a disassortative community structure with stronger connections between blocks than within blocks. If ϵ = 0, the three matrices are equal and the 9 networks have the same connection structure. Increasing ϵ differentiates the 3 sub-collections of networks. For δ-colSBM an δπ-colSBM , we add density parameters δ 1 = δ 4 = δ 7 = 1, δ 2 = δ 5 = δ 8 = 0.75 and δ 3 = δ 6 = δ 9 = 0.5.

α as = .3 +   ϵ -ϵ 2 -ϵ 2 -ϵ 2 ϵ -ϵ 2 -ϵ 2 -ϵ
Each of these configurations is simulated 30 times. We apply the strategy exposed in Section 6 and Appendix A of the main manuscript and evaluate the recovery of the simulated network partition.

Results

We assess the quality of our procedure by comparing the obtained partition of the collection of networks with the simulated one through the ARI index. As ϵ grows we are able to better differentiate the networks and do so almost perfectly on all colSBM setup. Note that Adding complexity slightly deteriorates the results as we recover the partition better for iid-colSBM and π-colSBM s than for δ-colSBM and δπ-colSBM .

S-4.4 Finding finer block structures

Finally, we perform a last simulation study in order to illustrate the fact that, for particular configurations, using a colSBM model on a collection of networks favors the transfer of information between networks and allows to find finer block structures on the networks. We consider the core-periphery structure configuration described in Equation (S-12) with ϵ = .4. In that case Q = 3.

We simulate a collection of 5 networks. 4 networks are of respective size (90,90,120,120). The last network is smaller with only 60 nodes and has a less marked structure (δ = .5) for the δ-colSBM and δπ-colSBM models.

Our goal is to recover the true connection structure of this last network X 5 . To do so, we compare the results obtained using either a standard single SBM on X 5 , or using the corresponding colSBM inferred with M = 2, 3 and 5 networks. We study Q in the various scenarii. In the simulation experiment, we obtained only Q = 2 or 3. The experiment is repeated 30 times. The results are depicted in Figure S-2.

For the 4 models of simulation, the simple SBM recovers 2 blocks most of the time. For iid-colSBM and π-colSBM , we always recover the 3 blocks while for the other case, we improve the ability to recover the true number of blocks when the quantity of information available from the other networks grows, either by augmenting the number of networks or by augmenting the number of nodes.

S-4.5 A note on model selection: dealing with networks of different sizes

We perform a simulation study to illustrate the fact that practitioner should be careful when dealing with networks of different sizes and should rely on BIC-L to state on the relevance of the common connectivity structure. We simulate 2 directed binary networks using sep-SBM with the following parameters: X as ∼ SBM 64 (3, π as , α as ) is drawn from an assortative community structure, while Table S-3: Average recovery of the simulated block memberships (ARI) for the assortative community networks (as) and the Erdős-Rényi networks (er) for each colSBM and the average difference in BIC-L between colSBM and the one used for the simulation (sep-SBM ).

X
to 640. The collection (X as , X er ) is simulated 20 times for each value of n er and inference is done for each colSBM . Our objective is to check if colSBM detects spurious structure on X er when n er is small and blurs the structure of X as when n er is large. We show the recovery of the grouping of the nodes into blocks for each network (ARI) and the difference in BIC-L with the best model (sep-SBM ) in Table S-3. In this setting, sep-SBM always detects the correct structure for all networks and no colSBM is ever selected by the BIC-L criterion. When n er is small, iid-colSBM , π-colSBM and δ-colSBM spuriously detect some structure on X er (low ARI), while δπ-colSBM is designed to correctly assign all the nodes of X er in one block. As n er gets larger, the absence of structure on X er is correctly detected by all models. When n er is much larger than n as , the structure of X as is blurred by the one X er while using iid-colSBM for the inference. The difference in BIC-L do not increase as fast for π-colSBM and δπ-colSBM compared to iid-colSBM and δπ-colSBM .

Figure 1 :

 1 Figure 1: Matricial view of 3 stream food webs. The species are reordered by blocks and blocks are ordered by expected out-degrees to emulate the trophic levels (bottom to top and right to left). The blocks have been obtained by fitting a SBM on each network separately. The shades of red depict the connectivity parameters α m .

  colSBM s with VEM starting from merged Ẑm sep-SBM (Q) (many initializations as a result of permutations within each Ẑm sep-SBM (Q)) -Keep the b fitted models with the best BIC-L for each Q while

Figure 2 :

 2 Figure 2: Estimated structure of the collection of 3 stream food webs with the four colSBM models. For each model, the matrix on the left is the estimated connectivity parameter α. The sizes of the blocks are proportional to the block proportion parameters π for iid-colSBM and δ-colSBM and proportional to the averages of block proportion parameters over the networks π m weighted by the number of species within the networks.. The barplot on the right depicts π (m) . The ordering is done by trophic level from bottom to top and right to left. For (δ-δπ)colSBM s we give δ below the barplot.

Figure 3 :

 3 Figure 3: Prediction of missing links and NA entries on stream food webs.

Figure 4 :

 4 Figure 4: Above: Clustering and connectivity structures of a collection of 67 predation networks from the Mangal database into 5 sub-collections. The length of the dendrogram is given by the difference in BIC-L to the best model. Below: Contingency table of the clustering found by π-colSBM and the different datasets from the Mangal database.

  Bernoulli, no explicit expression can be derived and one has to rely on a gradient ascent algorithm to update the parameters at each M-Step.

Figure S- 1 :

 1 Figure S-1: Partition of networks. ARI of the recovered partition of networks.

Figure S- 2 :

 2 Figure S-2: Finding finer block structures. Cumulative barplot of Q by the SBM s (blue) and the adequate colSBM (red) under the different simulation scenario. The number of blocks to be recovered is 3 and the darkest shade corresponds to Q = 3.

  as = (.4, .3, .3), α er = .25 , π er = (1).(S-13)

  er ∼ SBM ner (1, π er , α er ) is an Erdős-Rényi network with n er ranging from 10 er as er ∆ BIC-L as er ∆ BIC-L as er ∆ BIC-L as er ∆ BIC-L

		iid-colSBM	π-colSBM	δ-colSBM	δπ-colSBM
		ARI		ARI		ARI		ARI
	n 10 1	0	-5 1	.45	-7 1	0		-9 1 .9	-6
	20 1	0	-16 1	.35	-16 1	.15	-17 1 .75	-9
	40 .97 .1	-60 1	.45	-19 .99 .45	-41 1 .95	-9
	80 .94 .4	-95 1	.95	-14 1	.9	-45 1 1	-8
	160 .95 .9	-115 .99 1	-16 .99 1	-71 1 1	-12
	320 .9 1	-141 1	1	-24 1	1	-95 1 1	-18
	640 .7 1	-171 .99 1	-26 1	1	-132 1 1	-24

For the sake of clarity, we chose to use the terminology sub-collection, cluster and clustering for partitioning a set of networks while we use blocks or groups (and grouping) when referring to grouping nodes of a network into blocks.

Acknowledgments

The authors would like to thank Stéphane Robin for his helpful advice. This work was supported by a public grant as part of the Investissement d'avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH. This work was partially supported by the grant ANR-18-CE02-0010-01 of the French National Research Agency ANR (project EcoNet).

Model comparison

Estimation of Q and S

Parameter & Grouping accuracy (π-colSBM vs •) under π-colSBM under π-colSBM (mean ± sd) .

Results

For the model comparison, when ϵ α is small (ϵ α ∈ [0, .04]), the simulation model is close to the Erdős-Rényi network and it is very hard to find any structure beyond the one of a single block. As such, the iid-colSBM and π-colSBM models are equivalent and iid-colSBM is preferred to sep-SBM . We observe a transition when ϵ α = .08 where we become able to recover the true number of blocks Q = 4 and the support of the blocks given the true number of blocks. During this transition, the model selection criterion is about half of the time in favor of sep-SBM i.e. the model with no common connectivity structure between the networks.

From ϵ α = .16, we recover the true number of blocks and their support most of the time and the common structure obtained by the π-colSBM is found to be relevant. Note that when we are able to recover the true number of blocks, we are also able to recover their support almost every time.

For both the estimation of the parameters and the ARIs, the results mainly follow our ability to recover the true number of blocks, with the error of estimation of the parameters slowly decreasing from ϵ α = 0.12. ARI goes to 1 a bit faster than ARI 1,2 , denoting our ability to recover faster the real grouping of the nodes of each network than to match the blocks between the networks. This is directly linked with the detection of the true number of blocks and their support. Indeed, to get ARI 1,2 = 1, we need Rec( S, S) = 1 while the effective block number for each network is of only Q = 3, meaning that even with the wrong selected model we can still reach ARI = 1.

S-4.2 Capacity to distinguish π-colSBM from iid-colSBM

We aim to understand how well we are able to differentiate iid-colSBM from π-colSBM depending on the block proportions. To do so, we fix n m = 90 and Q = 3and set α and π as follows: with ϵ α = 0.16 and ϵ π taking 8 values equally spaced in [0, .28]. σ is a random permutation of the blocks. We simulate 30 different collections for each value of ϵ π .

Here again, we put in competition π-colSBM with iid-colSBM and sep-SBM and select a model if its BIC-L the greater than the two other ones. Then, for π-colSBM we compare Q to 3 and evaluate our ability to recover S. The results are provided in Table S-2. First notice that, since we chose ϵ π ≪ 1 3 , we do not simulate any empty block. As a consequence, the inference of the model is quite easy and we are able to recover the true number of blocks and the right support for the π-colSBM model almost always. When ϵ π = 0, π 1 = π 2 and the model reduces to iid-colSBM . This remark explains why iid-colSBM is preferred to π-colSBM when ϵ π < .2. As ϵ π increases, π-colSBM gets more and more selected, highlighting our capacity to recover the simulated structure.

S-4.3 Partitioning a collection of networks

The third simulation experiment aims to illustrate our capacity to perform a partition of a collection of networks based on their structure, as presented in Section 6 and Appendix A.

Simulation scenario For iid-colSBM , π-colSBM and δ-colSBM and δπ-colSBM , we simulate M = 9 undirected networks with 75 nodes and Q = 3 blocks. The block proportions are chosen as follows: .2, .3, .5) and for all m = 2, . . . , 9 π m = π 1 for iid-colSBM and δ-colSBM σ m (π 1 ) for π-colSBM and δπ-colSBM


Antoine Bonneau
email: antoine.bonneau@insa-lyon.fr

Frédéric Le Mouël

Fabien Mieyeville

Energy-efficient in-situ monitoring using on-device and distributed learning