{'fraction_non_alphanumeric': 0.04289536522914929, 'fraction_numerical': 0.015276619005149745, 'mean_word_length': 4.3321061512501915, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Block models for multipartite networks.Applications in ecology and ethnobiology

Modeling relations between individuals is a classical question in social sciences, ecology, etc. In order to uncover a latent structure in the data, a popular approach consists in clustering individuals according to the observed patterns of interactions. To do so, Stochastic block models (SBM) and Latent Block models (LBM) are standard tools for clustering the individuals with respect to their comportment in a unique network. However, when adopting an integrative point of view, individuals are not involved in a unique network but are part of several networks, resulting into a potentially complex multipartite network. In this contribution, we propose a stochastic block model able to handle multipartite networks, thus supplying a clustering of the individuals based on their connection behavior in more than one network. Our model is an extension of the latent block models (LBM) and stochastic block model (SBM). The parameters -- such as the marginal probabilities of assignment to blocks and the matrix of probabilities of connections between blocks -- are estimated through a variational Expectation-Maximization procedure. The numbers of blocks are chosen with the Integrated Completed Likelihood criterion, a penalized likelihood criterion. The pertinence of our methodology is illustrated on two datasets issued from ecology and ethnobiology.


Introduction
Networks have become fundamental tools in various fields, such as biology, ecological theory or sociology to name but a few. Statistical analysis aims to study the structure of these networks for instance by unraveling the clusters or communities of individuals (nodes) shaping the observed interactions (see Kolaczyk, 2009, for a review). The recent years have witnessed a growing interest for multi-layer networks which is a classical but wide terminology for a collection of related networks (see Pilosof et al., 2017, for a review of multi-layer networks in ecology for instance). Among others, multi-layer networks encompass dynamic networks, i.e. networks evolving with time (see Kim et al., 2018, for a review), multiplex networks when several types of relations are simultaneously studied on a common set of individuals (see Kéfi et al., 2016;Barbillon et al., 2016, for instance) and multipartite networks. Multipartite networks are a generalization of bipartite networks. In a bipartite network, the nodes (representing the interacting entities) are partitioned into two disjoint sets and an edge links a node from one set to a node from the other set (see Figure 1, left). In a multipartite network, the nodes are divided into more than two sets and edges link entities from different sets (see Figure  1, middle). In what follows, these pre-specified sets of nodes will be referred to as functional groups. Such multipartite networks arise in ecology when studying the interactions between several groups of species such as the interactions plant/pollinator, plants/ants, etc (Pocock et al., 2012;Dáttilo et al., 2016) or in biology when analyzing networks issued from multi-omics datasets involving proteins, etc. (see Pavlopoulos et al., 2018, for instance). Generalized multipartite networks are an extension of multipartite networks: the nodes are still partitioned into functional groups but the interactions may occur not only between different functional groups but also inside some of the functional groups. Ethnobiology (which is the scientific study of the relations between environment and people) provides such networks. One of the problematic of ethnobiology is to understand how social relations between individuals (seed circulation in our example) may structure and guarantee biodiversity in the cultivated crop species (see Thomas and Caillon, 2016, for example). As a consequence, the network of interest has two types of nodes -namely farmers and crop species-and we observe interactions between farmers and crop species (who grows what) and inside the group of farmers (seed exchange). Marketing also provides generalized multipartite networks when individuals are connected through social networks and also interact with goods through their on-line purchases.
Many statistical tools have been proposed to analyze simple or bipartite networks. However, when individuals or biological entities are involved in several networks, there is a strong need to propose statistical tools handling these networks jointly. These tools are an essential step towards a better understanding of the global interaction systems at stake. The aim of this paper is to propose A central issue in network analysis is to be able to cluster nodes sharing the same connectivity patterns. To do so, two approaches are possible, either relying on classical metrics detecting nestedness (Almeida-Neto et al., 2008), modularity (Barber, 2007), etc, or adopting a probabilistic mixture strategy. Since the original paper by Snijders and Nowicki (1997), probabilistic mixture models -classically referred as Stochastic Block models (SBM)-have proven their efficiency when aiming to cluster similar nodes based on their connectivity patterns without any a priori hypothesis about the patterns to be found (e.g. modularity, centrality, hierarchy). Latent Block Models (LBM) (Govaert and Nadif, 2003) are the extension of SBM for bipartite networks, resulting into a bi-clustering of the nodes.
In this paper, we propose a new probabilistic mixture model adapted to generalized multipartite networks. This model assumes that each functional group is partitioned into clusters gathering nodes sharing the same connection behavior on the basis of all the networks they are involved in. The clusterings are introduced through latent variables resulting into a mixture model on edges.
Inferring the parameters of such a model with latent variables is a complex task, since the likelihood can not be computed in a close form. We resort to a Variational Expectation-Maximization (VEM) algorithm to estimate the parameters and supplying the multi-clustering of the nodes. The crucial task of estimating the numbers of clusters in each functional group is tackled though a model selection procedure based on an adapted penalized likelihood criterion namely Integrated Classification Likelihood (ICL).
Related works Some extensions of standard descriptive tools -such as community detection-to more complex networks have been proposed in the literature (see for instance Yang and Leskovec, 2012;Gaskó et al., 2017). These metrics or descriptive tools have the drawback to only look for specific types of structures, contrary to clustering methods based on probabilistic mixture models. To the best of our knowledge, no such a probabilistic model has been proposed for generalized multipartite network. However, latent block models do exist for other kind of multi-layer networks. When the network is multiplex (i.e. when several types of interactions are studied on the same set of nodes), the SBM has been extended in Kéfi et al. (2016); Barbillon et al. (2016) by assuming a multivariate Bernoulli distribution to model the set of multiple edges between two individuals. Extending the SBM to a dynamic context where several snapshots of the networks at different time points are available, was proposed in numerous papers. A main difference between these papers is whether the evolution is considered as been discrete (Xu and Hero, 2014;Matias and Miele, 2017) or continuous (DuBois et al., 2013;Xin et al., 2017) in time.
Outline of the paper Section 2 is dedicated to the introduction of notations supplying a flexible tool to describe generalized multipartite networks. We also provide a description of the two datasets of interest and illustrate the notations on these specific cases. The block model for multipartite networks is described in Section 3. In Section 4, the variational inference the model selection procedure are presented including the practical implementation of the algorithm. Numerical illustrations of the robustness of the inference method are provided in Section 5. Finally, the statistical analyses of the two datasets with discussion are presented in Section 6. Perspectives are discussed in the last section. Note that this document is accompanied with a appendix sections divided into three parts: in the two first sections, the derivations of the update formulas for the VEM algorithm of the ICL criterion are provided; finally, detailed estimates on the two datasets are given.

Notations and data
A generalized multipartite network can be seen as a collection of networks involving Q functional groups: each network may be simple, i.e. describing the relations inside a given functional groups or bipartite, i.e. describing the relations between individuals belonging to two different functional groups. Let n q be the number of individuals in the q-th functional group (q = 1, . . . , Q). We index the collection of networks by pairs of functional groups (q, q ) (q and q in {1, . . . , Q}). The set E denotes the list of pairs of functional groups for which we observe an interaction network. For any (q, q ) ∈ E, the interaction network is encoded in a matrix X qq such that X qq ii = 0 if there is an edge from unit i of functional group q to unit i of functional group q , 0 otherwise. Each network may be binary (X qq ii ∈ {0, 1}) or valued (X qq ii ∈ N or R). X qq may be symmetric if the relation is non-oriented, non-symmetric otherwise. X = X qq (q,q )∈E encodes the generalized multipartite network. For each network, S qq is an additional notation which refers to the list of all the possible interactions.
• The dataset 2 comes from Thomas and Caillon (2016) and Thomas et al. (2015). They collected the oriented network of seed circulation between farmers -resulting in a non-symmetric adjacency matrix -and the crop species grown by the farmers, resulting in an incidence matrix. Noting q = 1 for the farmers and q = 2 for crop species we get E = {(1, 1), (1, 2)}. X 11 ii = 1 if farmer i gives seeds to farmer i (oriented relation), 0 otherwise and X 12 ij = 1 if farmer i cultivates crop species j, 0 otherwise.

A block model for generalized multipartite networks
In order to account for heterogeneity among individuals, we propose a mixture model explicitly describing the way edges connect nodes in the various networks. We assume that each functional group q is divided into K q blocks (or equivalently clusters). ∀q ∈ {1, . . . , Q} and ∀i ∈ {1, . . . , n q }, let Z q i be the latent random variable such that Z q i = k if individual i of functional group q belongs to cluster k. The random variables Z q i 's are assumed to be independent and such that: ∀k ∈ {1, . . . , K q }, ∀q ∈ {1, . . . , Q}, ∀i ∈ {1, . . . , n q }: with π q k ∈ [0, 1], ∀k = 1, . . . , K q and Kq k=1 π q k = 1, ∀q ∈ {1, . . . , Q}. We set Z q = (Z q i ) i∈{1,...,nq} , Z = (Z q ) q∈{1,...,Q} and π = (π q k ) k∈{1,...,Kq},q∈{1,...,Q} . Then, the nodes connect as follows: for any (q, q ) ∈ E, ∀(i, i ) ∈ S qq , F qq depends on the relation represented in X qq . We assume that for any (q, q ), F qq is either the Bernoulli distribution if the relation in X qq is binary, or the Poisson distribution if X qq ij is a counting or the Gaussian distribution if X qq encodes a continuous strenght of interaction. Equations (1) and (2) define the so-called Multipartite Block Model (MBM). Remark 1 Our model is a generalization of the SBM and the LBM. Indeed, Equations (1)-(2) reduce to the SBM if E = {(1, 1)} and to the LBM if E = {(1, 2)}. Our extension assumes that the latent structures Z are shared among the X qq i.e. if a functional group q is at stake in several X qq , the same Z q impacts the distributions of the corresponding interaction matrices. In other words, the clusters gather individuals sharing the same properties of connection in the collection of networks. Obviously, if each functional group appears in only one element of E, the MBM reduces to independent SBMs or LBMs. In terms of probabilistic dependence, conditionally on Z, the X qq ii are independent. However, Z being latent, their marginalization introduces dependence not only between the X qq ii but also between the matrices X qq (q,q )∈E . As a consequence, the clustering variables (Z q i ) q∈{1,...,Q},i ∈{1,...,nq} are dependent once conditioned by X.
For a given vector K = (K 1 , . . . , K Q ), let θ K = (α, π) be the unknown We have to perform two inference tasks: first, for a given vector K, estimating θ K and Z, second, selecting the right K.

Parameter inference and model selection
Let (X, Z; θ K ) denote the complete likelihood of the observations X and the latent variables Z for parameter θ K . Equations (1) and (2) lead to: The summation over Z in Equation (4) becomes quickly computationally intractable when K q 's and/or the n q 's increase. Moreover, we are interested in inferring the clusterings Z.

Variational EM algorithm
In such models with latent variables, the EM algorithm (Dempster et al., 1977) is a standard tool to maximise the likelihood, taking advantage of the simple form of log c (X, Z; θ K ). However, when conditioned by the observations X, the (Z q i ) are not independent (see Remark 1). As a consequence, the E step of the EM algorithm -which consists in the integration of log c (X, Z; θ K ) against P(Z|X; θ K )-has no explicit expression. In the context of the SBM and the LBM, the variational version of the EM (VEM) algorithm has proved to be a powerful tool for maximum likelihood inference (see Govaert and Nadif, 2008;Daudin et al., 2008). In the VEM, the problem of the E-step is tackled by replacing P(Z|X; θ K ) by an approximation R X,τ sought among distributions enforcing independence between the Z q i 's: Iterating steps (VE) and (M) leads to the maximisation of a lower boud of the likelihood: For the MBM, at iteration (t), the algorithm VEM is as follows: Step. Update of the parameters: Remark 2 Formula (5) on α qq kk is valid for the expectations of the Bernoulli, the Poisson and the Gaussian distributions. If F qq is a Gaussian distribution of variance v qq kk , then: The details of the VEM algorithm are provided in the Section A. As for any EMtype algorithm, its initialization is critical and will be discussed later. (θ K ,τ ) denotes the VEM estimates. A by-product of this algorithm is an approximation of the conditional distribution P(Z|X;θ K ) by R X,τ . Z is thus estimated by: The consistency of the VEM estimates has been established for the SBM by Bickel et al. (2013) and the LBM by Brault et al. (2017) while Mariadassou and Matias (2015) study the behavior of P(Z|X;θ K ) for the same models.

A penalized likelihood criterion
In practice, the number of clusters K = (K 1 , . . . , K Q ) is unknown and should be estimated. We adopt a model selection strategy where a model M K refers to the MBM with K clusters. Among many classical model selection criterion such as AIC, BIC and their variants, Biernacki et al. (2000) proposed the Integrated Classification Likelihood (ICL). ICL has proven its capacity to outline the clustering structure in networks in Daudin et al. (2008) (for simple networks), Keribin et al. (2014) (for bipartite networks) or Mariadassou et al. (2010) for valued networks. Its success comes from the fact that when traditional model selection criterion essentially involves a trade-off between goodness of fit and model complexity, ICL values not only goodness of fit but also clustering sharpness. Following the same line, we propose a modification of the ICL adapted to generalized multipartite networks: where andẐ has been defined in Equation (6). The better model is chosen asK = arg max K ICL(M K ). The penalization term pen(M K ) is made up of two parts: Q q=1 (K q − 1) log(n q ) corresponds to the clustering distribution and involves the number of nodes, while the other term depends on the size of α qq kk and the numbers of possible edges in each network |S qq |. The derivation of the ICL criterion is detailed in the Section B.

Practical algorithm
The practical choice of the model and the estimation of its parameters are computational intensive tasks. Assume that K q ∈ {1, . . . , K q }, then, ideally, we should compare Q q=1 K q models through the ICL criterion. For each model, the VEM algorithm has to be run starting from a large number of initialization points chosen carefully (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 VEM algorithm. The procedure we suggest is the following one. Given a current model , the m-th iteration of the procedure is written as follows.

Model selection strategy for MBM
• Split proposals. For any q such that K (m) q < K q , consider the model initializations by splitting any of the K · Compute the corresponding ICL(M (m+1)q + ) from formula (7). ).
• Set M (m+1) = arg max Initializing the VEM from the clusters obtained on a smaller or larger model is much more efficient than other strategies such as random initialization or spectral clustering. Note that the tasks at each iteration can be parallelized.
Simulations are implemented as vignette of the R-package available on Github (https://github.com/Demiperimetre/GREMLIN) in order to ensure reproducibility of the results.
Results for Scenario 1 Among the 100 simulated datasets, the true numbers of clusters are exactly recovered for 73 datasets as can be seen in Table 1. Three other estimated configurations are observed as detailed in Table 1, each of them corresponding tô K q = K q − 1. A scrutiny shows that cases where K 3 is estimated to 1 instead of 2 correspond to datasets where one simulated cluster (meaning the simulated Z 3 ) is reduced to 1 or 2 individuals.  We measure the ability of the procedure to recover the clusters by computing the Adjusted Rand Index (ARI by Hubert and Arabie, 1985) betweenẐ 1 and the simulated Z 1 . ARI compares two clusterings with a correction for chance. It is close to 0 when the two clusterings are independent and equals 1 when the clusterings are identical (up to label switching). The boxplots of the ARI's are plotted in Figure 2a. All the ARI are high (> 0.7), even whenK 1 = 6, which means that the clusters are always globally recovered and the structure of the multipartite network is well discovered. We also investigate the quality of the parameters estimation. For the 73 datasets such thatK = (7, 2, 2, 1), we compare the estimated parametersα to the simulation parameters (after a relabeling of the clusters to match the true clustering, if required). The results in terms of biais and Root Mean Square Error (RMSE) are reported in Table 2. We observe that the parameters are estimated without noticeable biases and with small RMSE. In particular, the null α qq kk 's are always recoverered.
Results for Scenario 2 The same analysis is done for Scenario 2.K = K = (3, 2) for 82 simulated datasets over 100. The unique alternative observed configuration isK = (2, 2). Here again, in order to assess the ability of the method to recover the clusters, we compute the ARI for Z 1 and Z 2 . These quantities are reported in Figure  2b for the two functional groups (left au right) and forK 1 = 3 (triangles) and K 2 = 2 (circles). As can be noticed, the clusterings are well recovered when the number of clusters are exactly recovered (triangle), with a noticeable number of cases where the clusters are exactly recovered (ARI= 1). As expected, we observe a lower ARI whenK 1 = 2 (circle), but with globally high values.
We also pay attention to the estimation of the parameters of connection α qq kk . Biais and RMSE are reported in Table 3. Once again, the biais and RMSE are small.

Ecology: interactions plants / animals
Dataset The dataset -compiled and conducted by Dáttilo et al. (2016) at Centro de Investigaciones Costeras La Mancha (CICOLMA), located in Mexico-involves three general types of plant-animal mutualistic interaction: pollination, seed dispersal by frugivorous birds, and protective mutualisms between ants and plants with extrafloral nectaries. The dataset -which is one of the largest compiled so far with respect to species richness, number of interactions and sampling effort-includes n 1 = 141 plant species, n 2 = 173 pollinator species, n 3 = 46 frugivorous bird species and n 4 = 30 ant species, inducing 753 observed interactions of which 55% are plant-pollinator interactions, 17% are plant-birds interactions and 28% are plant-ant interactions.
Inference We run the procedure described in Subsection 4.3 starting from several automatically chosen initial points K (0) , with numbers of clusters bounded between 1 and 10. The computational time on an Intel Xeon(R) CPU E5-1650 v3 3.50GHz x12 using 6 cores is less than 10 minutes. The ICL criterion selects 7 clusters of plants, 2 clusters of pollinators, 1 cluster of birds and 2 clusters of ants. The estimated parameters are reported in Tables 5 and 6 in the supplementary material. They are really similar to the ones used for the simulation study and provided in Equation 8. The resulting mesoscopic view of the multipartite network is plotted in Figure 3. Discussion From Figure 3, we conclude that plants of Clusters 7 and 2 only interact with ants (plants of Cluster 7 attract more ant species belonging to Cluster 1). The plants of Clusters 3 and 6 are only in interaction with birds, the difference between the two clusters being due to the strength of the connection. The difference between the two clusters of pollinators derives from the existence of Cluster 1 of plants.
In order to illustrate the contribution of our method, we also analyze each bipartite network separately (using an LBM) and compare the results in terms of clustering. The clusterings are compared through the Adjusted Rand Index (ARI): if ARI = 1 then the clusterings are equal (up to a label switching transformation). The ARIs are given in Table 4.
Using standard LBM, we obtain 2 clusters of ants, 1 cluster of birds and 3 clusters of pollinators. The clusterings of ants and birds are not modified by the ecosystemic approach, their ARI being equal to 1. The clustering of the pollinators is slightly modified, going from 3 clusters to 2 clusters but the additional block only contains few individuals, thus leading to a high ARI. Since the plants functional group is involved in the three bipartite networks, we obtain 3 clusterings when analyzing them separately. These three clusterings are -as expected-very different from our MBM clustering (the ARIs being respectively equal to 0.118, 0.415 and 0.163, see Table 4). When aiming at proposing a clustering taking into account the 3 bipartite networks, one may adopt a naive strategy by combining (by intersection) these three clusterings. We then obtain 12 blocks of plants and the ARI with the MBM clustering is 0.617. However, this number of clusters (12) is too large with respect to the model selection criterion. Our MBM clustering is a parsimonious trade-off.
Finally, Figure 3 highlights the central role played by the pollinators in the complete network also involving ants and seed dispersal birds. Our ecosystemic approach unables us to unravel such a central position of pollinators in general with respect to ants and seed dispersal birds, which would have been impossible when performing separate analysis.

Seed circulation and crop species inventory
Dataset Seed circulation among farmers is a key process that shapes crop diversity Pautasso et al., 2013). Data on seed circulation within a community of first-generation migrants (30 farmers) were collected by a field survey in the island of Vanua Lava in the South Pacific archipelago nation of Vanuatu. A farmer is considered as a giver for another farmer if he/she has given at least one crop since they arrived in the new settlement site in Vanua Lava. It results in a connected and directed network of seed circulation. Besides the circulation network, inventory data for each farmer were also collected. They consist in the list of crop landraces they grow. This was aggregated at the species level, leading to 37 crop species. These inventory data can be seen as a bipartite network. The seed circulation data were analyzed in Thomas and Caillon (2016) and the inventory data were analyzed in a meta-analysis in Thomas et al. (2015). On the basis on the joint modeling we propose in this paper, we aim to provide a clustering on farmers and crop species on the basis of the seed circulation network and the inventory bipartite network.
Inference Three clusters of farmers and two clusters of crop species were selected. The inferred parameters are given in Table 7 and a a mesoscopic view is displayed in Figure 4. The clusters were renumbered to make them correspond to the probability of connection: the larger cluster number, the larger marginal probability of connection for the members of the cluster.
Discussion The discovered clusters are straightforwardly interpretable: Cluster 3 gathers farmers who circulate seeds within the cluster and give to the two other clusters, Cluster 2 circulates seeds within the cluster contrary to Cluster 1 who only receives from Cluster 3; the two clusters of crop species are Cluster 2 with more common crop species and Cluster 1 with other species. Clusters 3 and 2 of farmers grow crop species from Clusters 1 and 2 whereas Cluster 1 of farmers grows mainly crop species from Cluster 2. It turns out that Cluster 3 gathers mainly the first migrants and Cluster 1 the last migrants. The pattern of connection is then explained by the fact that first migrants helped the others to settle by providing seeds. Moreover, the first migrants had more time to collect more crop species to grow. In order to compare the clusterings obtained by the MBM and the ones obtained from the circulation network (clustering on farmers from SBM) and the inventory network (clusterings of farmers and plants with a LBM), we compute the ARI. The clusterings on crop species remain quite close since their ARI is equal to 0.891. However, the clusterings on farmers are quite different (ARI smaller than 0.3), indeed the MBM shall make a trade off between the circulation and the inventory for farmers. To ease the comparison between clusterings on farmers, the same renumbering rule was applied for all the different clusterings so that the larger cluster numbers correspond to larger marginal probability of connection. Figure 5 is an alluvial plot which compares the three obtained clusterings of farmers. It shows how the trade-off is made between the two stand-alone clusterings in the MBM clustering. It appears quite obvious that Cluster 1 given by the MBM gathers only farmers from Cluster 1 in the seed circulation network and from Clusters 1 and 2 from the inventory data since this cluster aggregates farmers with fewer connections and who grow less crop species than the others. The same kind of observation can be made for Cluster 3 given by the MBM which aggregates farmers who were in the cluster of the most connected farmers and in the two clusters of farmers who grow more seed.

Discussion
In this paper, we proposed an extension of LBM and SBM which can handle multipartite networks, resulting in the so-called MBM. Multipartite networks encompass a lot of various situations such as the two examples dealt with in the paper. Besides, MBM can be also useful for many other contexts with different multipartite structures. The main limiting factor of the parameter inference and model selection methods we propose, is the size of networks. Indeed, the inference algorithm is suitable for networks up to 1000 nodes in order to keep computational time reasonable. If willing at studying larger networks, one should develop adapted inference algorithms. Several extensions can be thought of with no additional significant difficulty. For instance, covariates can be taken into account by writing P(X qq ii = 1|Z q i = k, Z q i = k ) = φ α qq kk + y qq ii β where y qq ii are the covariates describing the pair of individuals (i, i ). A more challenging extension could be to include extra parameters in the model to account for degree heterogeneity in the vein of the degree-corrected SBM (Karrer and Newman, 2011). However, this should lead to cluster individuals independently of their degrees which is not always desirable. For instance in the seed circulation example, the in-degrees and Figure 5: Alluvial plot comparing the clustering on farmers obtained from an SBM on the circulation network (clExc), an LBM on the inventory network (clInv) and the MBM (clMBM) on both networks. The cluster numbers are related with the probability of connection, the larger cluster number, the larger marginal probability of connection (between farmers for clExc, between farmers and crop species for clInv). out-degrees of farmers are key factors to explain their social structure. In the interaction networks in ecology, distinguishing generalist (high degree) from specialist (low degree) species is a part of the expected outcomes of the clusterings. More recent works (Zhu et al., 2014) propose to use node degrees and edge orientations to help the clustering while tolerating heavy-tailed distribution of degrees. Incorporating these variants may be of interest in the MBM framework and could lead to further developments. In a more general perspective, the study of ecological or sociological interactions supplies a wide variety of complex networks such as multilevel networks or multi-layer networks (Pilosof et al., 2017;Lazega and Snijders, 2015). Some of them can be treated as multipartite networks and then by a MBM (possibly by incorporating some extensions discussed above). The others require the development of suited models which might also rely on a latent variable modeling. They will be the subject of future works.

A Variational EM for MBM: details
The variational version of the EM algorithm maximizes a lower bound of the observed likelihood. More precisely, let R τ ,X be any probability distribution on Z, we define I θ (R τ ,X ) as: where KL is the Kullback-Leibler divergence and H(R τ ,X ) is the entropy of R τ ,X . The inequality in (11) derives from the positivity of the KL divergence. The equality I θ (R τ ,X ) = log (X; θ) holds iff R τ ,X (Z) = P(Z|X; θ) which results in the classical identity of the EM algorithm (Dempster et al., 1977). The principle of the VEM algorithm is to choose R τ ,X in a family of distributions P parametrized by τ , such that the conditional expectation E R τ ,X [log c (X, Z; θ)] can be computed explicitly. Iteration (t) of VEM consists in the two following steps: Step The variational EM generates a sequence (θ (t) , τ (t) ) t≥0 increasing the lower bound I θ (R τ ,X ) of the likelihood log (X; θ).

Choice of P
The key point of the procedure is the choice of P, making the calculus tractable but rich enough to be a good approximation of the true conditional distribution P(Z|X; θ). Following Govaert and Nadif (2008) and Daudin et al. (2008), we adopt the mean-field strategy (Jaakkola, 2000) and chose P as: where h Kq (·; ξ) is the density of a 1 trial -multinomial distribution of parameter ξ ∈ T Kq , i.e. R X,τ is such that the latent variables Z are independent and From this particular choice of P, we derive the VE-step and M-step for the MBM in what follows.
Explicit expression for I θ (R τ ,X ) Using the expression of the complete loglikelihood given in Equation (3), we obtain: for weigthed Poisson networks. log φ(X qq ii , µ qq kk , σ qq kk ) for Gaussian networks.
has to be carefully calculated, when i = i . To that purpose, let us introduce the following notations : • ∀q, E q = {q ∈ {1, . . . , Q} | q = q and (q, q ) ∈ E} . E q is the set of incidence matrices involving the functional group q.

Clustering and initializations
We denote by θ and τ the resulting estimates. The estimated clustering is the maximum a posteriori (MAP) estimate: ∀q ∈ {1, . . . , Q}, ∀i ∈ {1, . . . , n q }, B Derivation of the ICL criterion As exposed in Section 4.2 of the paper, we resort to the ICL criterion to perform model selection. The ICL is an asymptotic approximation of the integrated marginal complete likelihood. We supply here the details of the calculations.
Explicit expression of the marginal complete likelihood The ICL being a based on a Bayesian model selection criterion, we first set the following prior distribution: α qq kk ∼ B(a, a) and (π q 1 , . . . π q Kq ) ∼ Dir(b, . . . , b).
where B(a, a) denotes the Beta distribution and Dir(b, . . . , b) is the Dirichlet distribution. By definition, log m c (X, Z; M K ) = log c (X, Z; θ K )π(θ K ; M K )dθ K .
Taking advantage of the conditional independences in the model defined by Equations (1) and (2) combined with the independence of the parameters in the prior distribution, we can decompose log m c into the following sum: log m c (X, Z; M K ) = (q,q )∈E log f (X qq |Z q , Z q ; (α qq ))π(α qq )dα qq + Q q=1 log f (Z q ; π q )π(π q )dπ q .
Using the fact that f (Z q ; π q ) = Kq k=1 (π q k ) N q k with and the conjugacy of the Dirichlet prior distribution, we easily deduce the following formula: f (Z q ; π q )π(π q )dπ q = Γ(bK q ) Γ(b) Kq Kq k=1 Γ(N q k + b) Γ(n q + bK q ) where Γ is the Gamma function. Now, we can reformulate f (X qq |Z q , Z q ; α qq ) as: With the beta prior distribution on each α qq kk , we get: f (X qq |Z q , Z q ; (α qq ))π(α qq )dα qq = where N q k has been defined in Equation (18) and S qq kk and N qq kk in Equation (19) Asymptotic approximation Using the same arguments as in Daudin et al. (2008) and Brault (2014), we obtain the following asymptotic approximation. Assume that ∀q ∈ {1, . . . , Q}, n q → ∞, then : The first term comes from the application of the Stirling formula to the Gamma function when approximating f (Z q ; π q ). The second term comes from a BIC approximation of the part f (X qq |Z q , Z q ; α qq ). Obviously, the parameters of the prior distribution (a, b) disappear from the formula since asymptotically the importance of the prior distribution vanishes.   Figure 3. π 1 1 0.4675π 2 1 0.06 2 0.1606 2 0.94 3 0.1351π 3 1 0.1 4 0.0784 * 2 0.9 5 0.1061π 4 1 1 6 0.0142 7 0.0380 Table 6:π q k for MBM on the plant/animals interactions dataset Table 7: Estimated parameters for MBM on the seed circulation dataset:π q k are in the first row and column, other rows and columns contain the estimates α qq kk . α qq kk identified by * are not plotted in Figure 4.