Generative model for reciprocity and community detection in networks

We present a probabilistic generative model and efficient algorithm to model reciprocity in directed networks. Unlike other methods that address this problem such as exponential random graphs, it assigns latent variables as community memberships to nodes and a reciprocity parameter to the whole network rather than fitting order statistics. It formalizes the assumption that a directed interaction is more likely to occur if an individual has already observed an interaction towards her. It provides a natural framework for relaxing the common assumption in network generative models of conditional independence between edges, and it can be used to perform inference tasks such as predicting the existence of an edge given the observation of an edge in the reverse direction. Inference is performed using an efficient expectation-maximization algorithm that exploits the sparsity of the network, leading to an efficient and scalable implementation. We illustrate these findings by analyzing synthetic and real data, including social networks, academic citations and the Erasmus student exchange program. Our method outperforms others in both predicting edges and generating networks that reflect the reciprocity values observed in real data, while at the same time inferring an underlying community structure. We provide an open-source implementation of the code online.


I. Introduction
Reciprocity in directed networks, i.e., the tendency of a pair of nodes to form mutual connections between each other [1], is an important feature of many social relationships. Its impact ranges from affecting the development of exchange and power to determining the emergence of trust and solidarity [2,3]. Behavior of this kind has also been found in many kinds of networks that reflect human and institutional interaction, e.g., the world wide web, online dating, interfirm contracts, journal citations and email communication [4][5][6][7][8].
Among the various network modeling approaches, that of probabilistic generative models enable us for a rigorous theoretical foundation within the framework of statistical inference, as well as a flexible incorporation of domain knowledge in the modeling assumptions. Here, we consider a latent variable model, a probabilistic approach that contains latent and observed variables. The latent variables encode hidden patterns in the data, such as community memberships, and determine the probability of ties between nodes. For instance, knowing which communities two nodes belong to helps determine the likelihood of their interaction.
While in some simple cases, community structure may explain the tendency toward reciprocation [9], this mechanism may not be enough to capture more complex scenarios. Indeed, many generative models with community structure fail to reproduce the values of reciprocity observed in real networks, as we discuss in more details later. Conversely, several models aimed at capturing reciprocity do not account for community structure [10,11]. It is reasonable to expect that the mechanism regulating the existence of interactions can be influenced by both patterns of communities and reciprocity. In addition, communities are often interpretable objects and may correspond to functional unit, hence the value of including them in the model formulation. Incorporating reciprocity as well as community structure into a coherent latent variable model comes with the main challenge of relaxing the conditional independence assumption between edges, a common assumption in generative models to ease mathematical derivations. In addition, this task requires properly capturing conditional probabilities, as we describe later. Inspired by these insights, we propose a novel probabilistic latent variable approach to model networks that preserves the benefits of generative models, while capturing both community structure and reciprocity.
Models for reciprocity and latent community structure have largely been developed independently of one another, and only a handful of works have hinted at incorporating them into a unique framework. For instance, Garlaschelli and Loffredo [12] point towards a possible relationship between their model for reciprocity and general hidden variable models. Most notably, the pair-dependent stochastic block model of Holland et al. [9], well explained also by Wasserman and Anderson [13], holds assumptions similar to ours, in that it models jointly pairs of edges, which they call dyad vectors. While a seminal work, it is, nevertheless, limited to hard membership and unweighted networks; hence the likelihood function that they propose substantially differs from the likelihood represented by our model. One practical aspect of our choice for the likelihood is that parameters' inference in our model is optimized to fully exploit the sparsity of the dataset and is scalable to large network sizes.
Reciprocity is often modeled by means of exponential random graphs [10,11,14,15], where it is treated as a measured network property that needs to be reproduced (often together with other network properties like the degree) by sampling networks using statistical mechanics principles, e.g., maximum entropy. The approach presented in this work significantly differs from the previous studies in that we include latent variables, such as community membership, as a mechanism to determine edge formation. However, in the case of exponential random graphs, possible group structures are not given a priori as the latent parameters; instead, they can only be estimated a posteriori on the sampled networks. More broadly, our approach is that of generative models, which incorporate a priori community structure by means of latent variables, and these are inferred from the observed interactions [16,17]. However, in these generative models reciprocity is not explicitly included as a mechanism for tie formation, thus these models often fail to reproduce the observed reciprocity values of real networks. Consequently, a generative method whose latent variables describe both reciprocity and community memberships is needed.

II. Relaxing the conditional independence assumption
A possible explanation for the practical deficiency of generative models with communities to reproduce observed reciprocity values is the common assumption of conditional independence between edges, which makes the problem both analytically and computationally more tractable. This assumption states that the likelihood of a directed tie between two nodes depends only on their community membership (and other possible model parameters), but not on the existence of the reciprocated edge. This might be too strict of an assumption to capture the feature of reciprocity, where it is reasonable to expect that the existence of an edge in one direction should also be conditioned on the existence of an edge in the opposite direction. For instance, if an author i has cited another author j, this might predict the probability of j also citing i. At the same time, knowing the communities that the authors belong to, could also help estimating this probability. Mathematically, this can be translated to relaxing the assumption of conditional independence, which is the approach we take here.
Formally, we represent interactions between N individuals as a weighted asymmetric matrix A, with entries A ij being the number (or weight) of interactions from i to j; for instance, the number of favors or services that i does for j, or the number of times that i has endorsed j, e.g., as paper citations. Our model assigns a joint likelihood P (A ij , A ji |Θ) to edges involving the same pairs of nodes (i, j), given some set of latent parameters Θ. Specifically, we assume the likelihood of a network to factorize as: This is fundamentally different from the prevalent approaches in generative models, where, typically, one assumes that individual edges are conditionally independent given the network parameters, i.e., P (A|Θ) = i,j P (A ij |Θ).
Notice that edges involving different pairs of nodes remain conditionally independent as in standard approaches. Equivalently, in terms of the conditional distribution of an individual edge P (A ij |A ji , Θ), we assume that this can be different than its marginal distribution P (A ij |Θ). To the extent of our knowledge, this assumption has never been deeply questioned, except for a few works [18,19]. As firstly pointed out by Hoff [20], there are theoretical groundings for this assumption to hold in common scenarios, due to generalizations of de Finetti's theorem by Aldous [21] and Hoover [22] (see [19] for a detailed discussion). They show that, for exchangeable graphs, i.e., in networks without any natural order between nodes (which is often the case), the joint probability function of the adjacency entries can be properly represented using latent variables on nodes and pairs. In other words, the joint can be factorized as a product on edges, given the latent variables.
However, in the case of directed networks, where the adjacency matrix is asymmetric, as in our case, a precise representation can only be obtained using Eq. (1). While standard conditionally independent models can in principle arbitrarily well approximate the whole network distribution [23], in practice it is not known how stateof-the-art models perform on this regard. To effectively model reciprocity, we relax the assumption of conditional independence and include the pairwise dependencies of two directed edges between pairs of nodes; such minimal relaxation is required to effectively model reciprocity. We compare results against standard conditionally independent models in terms of various performance metrics on both synthetic and real data.

III. The community-reciprocity model
To fully specify the joint likelihood in Eq. (1), we need to characterize conditional distributions and onepoint marginals like the distribution P (A ij |A ji , Θ) and P (A ij |Θ). Here, we aim at capturing reciprocity, hence we assume that observed interactions exist because of two types of contributions: i) the communities that nodes belong to, as in general community detection frameworks like the stochastic block model [9], and ii) the fact that an individual that receives a directed interaction is more likely to reciprocate. In order to construct a model flexible enough to capture weighted networks and overlapping communities, we utilize a mixed-membership approach, similar to [16,17], to model how communities regulate edge formation.
Given the adjacency matrix A, our goal is to find community memberships of nodes and the magnitude of the reciprocity effect in the network, i.e., Θ. Bringing the contributions of reciprocity and community structure to-gether, we model the conditional probability of A ij given A ji as drawn from a Poisson distribution with mean We denote with Θ = (u, v, w, η) the latent parameters that we want to infer. The parameters u ik , v ik are entries of K-dimensional vectors u i and v i , the out-going and in-coming communities respectively; w kq are the entries of a K × K affinity matrix, which regulates the structure of communities, e.g., assortative when its diagonal entries are greater than off-diagonal entries, in this case edges are more likely between nodes in the same community; η is the reciprocity parameter, and it regulates the impact of observing A ji to predict the existence of A ij . We omit from it the number of communities K, as in this work we assume this as given.
When unknown, as in our experiments with real data, we estimate it by using cross-validation.
Notice that λ ij includes separate contributions from both community parameters and reciprocity coefficient. It assumes additive contributions: we can have zero contribution from one term and still observe the existence of an edge because of the other term. If both are non-zero, their total impact sums up. This is conceptually different than a multiplicative contribution, a possible modeling choice that we do not explore here. Intuitively, an edge with weight A ij exists if i and j belong to compatible communities (compatibility is regulated by the affinity matrix) or because of the reciprocity effect of observing the opposite edge A ji . For instance, an author might cite another one because they belong to the same community (e.g., a research sub-field) or because she was cited by the other on a previous publication.
Finally, as we need positive λ ij , we assume η ≥ 0. This restricts the model to have positive reciprocity contribution, i.e., receiving an in-coming edge can only boost the likelihood of the corresponding out-going edge, but not decrease it. Although this assumption could be limiting in certain contexts, it nevertheless applies to several relevant scenarios, in particular to the cases we study here. Relaxing this assumption, and suitably modifying the underlying theoretical model, is left for future works.
Our model specifies conditional probabilities, however, we do not assume the existence of a consistent joint distribution. In fact, finding a closed-form for the joint in Eq. (1), consistent with our proposed conditional, requires specifying a marginal probability function and then enforce consistency equations like Aji P (A ij |A ji , Θ) P (A ji |Θ) = P (A ij |Θ). Depending on the choice of this marginal, enforcing consistency might be non-trivial, as it may require performing intractable marginalization. Early formalizations of the consistency between conditional and joint distribution has been provided, in a seminal work, by Besag's Auto-Poisson models [24]. In the context of graphical models, a few models specify conditional Poisson distributions [25,26], but without considering latent variables. In the absence of a closed-form joint distribution, we adopt a tractable pseudo-likelihood approach [24], where instead of optimizing the exact likelihood of Eq. (1), we consider the approximation: (4) which is available in closed-form as it requires only the conditional probabilities, which we specified above. The equality holds only when A ij and A ji are conditionally independent, the common assumption in network generative models, as in that case P (A ij |A ji , Θ) = P (A ij |Θ). This is not our case since we relax this assumption, and Eq. (4) is only an approximation. This approach has also been considered in dyadic-dependent models [27], for community detection in networks [28], and for local Poisson graphical models [25]. A visual overview of our model is shown in Fig. 1.

IV. Inference with expectation-maximization
The goal is to find the community and reciprocity parameters, i.e., Θ, given the adjacency matrix. Defining L ps ij (Θ, A ji ) = log P (A ij |A ji , Θ) and neglecting the factorial term which is independent of these parameters, we have the log-pseudo-likelihood: We aim at maximizing this quantity, but the presence of the logarithmic term makes this maximization difficult. However, using a variational approach by means of Jensen's inequality, it can be shown (see Appendix D 1) that maximizing L ps (Θ) is equivalent to maximizing with respect to Θ, ρ = ρ (1) , ρ (2) , and φ, where are the variational distributions over the parameters.
Constraints on the parameters can be arbitrarily added, e.g., k u ik = k v ik = 1, by incorporating Lagrange multipliers inside Eq. (5), and repeating similar calculations. In our numerical experiments, we consider both constrained and unconstrained cases.
We can perform this optimization by alternatively updating the various parameters, with an expectationmaximization (EM) algorithm. At each step, one updates ρ and φ using Eqs. (7)-(8) (E-step) and then maximizes L ps (Θ, ρ, φ) with respect to Θ by setting partial derivatives to zero (M-step). This iteration is repeated until L ps convergences. The whole routine is described in Algorithm 1 and the detailed derivations are in the Appendix D. This algorithm is computationally efficient and scalable to large system sizes as it exploits the sparsity of the dataset. Indeed, all the updates involve in the numerator sums over A ij , hence only the non-zero entries count, giving an algorithmic complexity of O(M K 2 ).

Algorithm 1 CRep: EM algorithm
network-affinity matrix w = [w kq ]; reciprocity parameter η. Initialize u, v, w, η at random. Repeat until L ps convergences: 1. Calculate ρ (1) and φ (E-step): i) for each node i and community k update memberships: ii) for each pair (k, q) update affinity matrix: iii) update reciprocity parameter: Note: γ u i , γ v i are quantities that are defined differently for constrained and unconstrained values of u i and v i . In the constrained case, they correspond to Lagrange multipliers; see Appendix D 2.

V. A benchmark generative model with communities and reciprocity
So far we have focused on recovering the model parameters given the data, i.e., the inference. In this section, instead, we propose a benchmark probabilistic generative model to generate synthetic data with intrinsic community structure, and a given reciprocity value. It takes as input a set of membership vectors, u i and v i , affinity matrix w, and reciprocity parameter η; the output is a directed network with adjacency matrix A. In this formulation, edges between a given pair of nodes are generated stochastically; one edge being generated first and independent from the other, while the formation of the opposite edge depends on how the first was drawn. The pairs of edges are conditionally independent from each other. Formally, we aim at sampling pairs of edges from Eq. (1), which can be done by selecting a marginal P (A ij |Θ) and a conditional distribution P (A ji |A ij , Θ). By assuming a Poisson conditional as in Eq. (2) and a Poisson marginal, our model would reduce to a standard (conditionally independent) generative model with communities in the case of zero reciprocity parameter. Even though with this choice the joint is computationally intractable, this is not an issue, as we do not aim to use the joint to compute quantities analytically, but rather focus on sampling from it. Formally, given the input set of latent variables Θ = (u, v, w, η), we draw a pair (A ij , A ji ) consistently with the joint P (A ij , A ji |Θ), in a two-step sampling routine: 1. Select with a coin-flip one direction, (i, j) or (j, i).
Say we select (i, j).

Sample
where is the mean of the marginal distribution such that it is consistent with the joint and the con- using the previously extracted value of A ij .
The Poisson distribution may generate multiple edges between a pair of nodes, so this model may create multigraphs. This is consistent with the interpretation that A ij is the number, or total weight, of links from i to j. If we wish to generate binary networks where A ij ∈ {0, 1}, we use the fact that the Poisson and Bernoulli distributions become close in the sparse limit. To enforce sparsity, it is sufficient to multiply the λ 0 ij by a constant ζ, as the m ij in Eq. (10) will also be automatically rescaled by the same quantity. The constant can be fixed by choosing a value for the expected number of (weighted) edges: Imagine now a practitioner willing to control for the relative contribution of community and reciprocity in generating edges. Our model naturally allows this possibility, as this tuning is encoded by η. To see this explicitly, we calculate the fraction of edges generated by community effects only and introduce the cr ratio variable as following: where we used Eq. (10) to rewrite the denominator. Thus, by varying η in the input, one automatically tunes the interplay community vs reciprocity: η close to 0 gives a network whose edges depend mostly on the community structure imposed by the membership vectors; instead, η close to 1 results in a network with lower impact of community structure, i.e., reciprocity has also significant impact on the edge formation. Notice that it is not possible to have a contribution purely due to reciprocity, as this phenomenon implicitly requires the existence of another mechanism to produce one of the two possible edges, here the community structure. This can also be seen by observing that Eq. (10) can be rewritten as m ij = λ 0 ij + η 1−η 2 η λ 0 ij + λ 0 ji ; while the first term only depends on communities, the second term depends on both communities and reciprocity and they cannot be separated independently.
Having presented how our model can be used to generate synthetic data, we now proceed in describing how our model relates to observable network properties and how it can be used to predict reciprocated edges.

VI. Predicting network reciprocity
In directed networks, reciprocity r is usually defined as the fraction of edges that are reciprocated [1], although other definitions exist to capture this feature [15,29]. With our probabilistic model, we can compute the expected value of a related quantity which corresponds to reciprocity in the case of binary adjacency matrices. A natural question is thus how this observable quantity is related to the reciprocity parameter η. In fact, we show that, provided some assumptions for the second moment E A 2 ij and considering an approximation with Taylor expansion (see Appendix D 4), η is a lower bound for it: The tightness of this bound depends on the latent variables through λ 0 ij , (implicitly) m ij , and m ji . Empirically, we find that in the majority of the experiments the bound is very tight, i.e., E [r w ] ≈ η and the other terms in Eq. (16) are much smaller than η in models with the conditional independence assumption, such as our proposed model with η = 0, where E [r w ] = i,j mij mji i,j mij . In fact, in these models, the term i,j m ij m ji is often very small -we show empirical evidence of this later. Therefore, even in networks with high reciprocity, models with conditional independence assumption could poorly reproduce the term. This empirical result also seems to indicate that the pseudo-likelihood approximation of Eq. (4) is relatively good in our datasets. The practical indication for practitioners is that networks generated by models with the conditional independence assumption have reciprocity values significantly different from those observed in real data.

VII. Predicting reciprocated edges
The dependence structure between pairs of edges should allow us to predict the existence of a reciprocated tie if an edge in the opposite direction is observed, such as the citation of a paper if an author has been cited before by someone else. This is a kind of link prediction task, which lets us test the dependence assumption. It is also a principled way of comparing the accuracy of various generative models for any real network where no ground truth for the latent variables is known [30].
Conditional edge prediction can be formulated as follows: what is the probability of an edge i → j conditioned on observing the opposite existing edge (or non-existing edge) j → i ? Our model naturally outputs this conditional probability. In contrast, a generative model that assumes conditional independence between edges is not capable of exploiting this extra information. It could only estimate marginal probabilities that do not depend on observing the opposite edge as it uses only the parameters such as community memberships and the affinity matrix. Our model is not capable of fully estimating marginal distributions but nevertheless can estimate its expected value as in Eq. (10). This is often the main quantity used in prediction tasks, as it plays the role of a score for estimating the entries A ij . Therefore, with our model we can also predict regular edge existence, where we simply aim at predicting an edge without any extra information but the inferred parameters.
In our experiments below, we test various generative models for both regular and conditional edge prediction by using 5-fold cross-validation. Specifically, we divide the dataset into five equal-size groups (folds) and give the models access to four groups (training data) for learning the parameters; this contains 80% of the possible pairs of nodes in the network. One then predicts the existence of edges in the held-out group (test set). As performance metrics, we measure the AUC on the test data, i.e., the probability that a randomly selected edge has higher expected value than a randomly selected non-existing edge. We compute both the regular AUC, by using as score the expected value E P (Aij |Θ) [A ij ], and the conditional AUC (AUC−cond), which uses E P (Aij |Aji,Θ) [A ij ] as the score, i.e., the expected value over the conditional distribution. The latter can only be computed for our algorithm, as for the others the marginal distribution is the same as the conditional, and thus the two AUC values coincide, see Appendix B for more details.

A. Results on real and synthetic data
We now demonstrate our model by applying it to both real and synthetic data. In the real-world datasets available to us, we only have a directed network of observed interactions, i.e., there is no available ground truth for the actual membership and reciprocity parameters. Consequently, their relative contributions in edge formation cannot be tuned. Thus, we first validate our model and competing algorithms on synthetic data produced with different generative models. We test the ability of these models to: i) generate sample networks that replicate relevant network quantities such as reciprocity, similar to the observed values on the input networks; ii) perform edge prediction tasks. We then investigate our model's performance on real-world datasets.
In the tests below, we use our model in various ways: the constrained version with constraints on the membership parameters u and v such that k u ik = k v ik = 1, ∀i (CRep), the non constrained version (CRepnc), and our model with η = 0 (CRep0), i.e., without considering the reciprocity effect. For comparison, we use two generative models with latent variables: a community detection-only generative model with a Maximum Likelihood approach [16] (MT), which was the inspiration for the building block of our model in the case η = 0, and a Bayesian Poisson matrix factorization (BPMF) commonly used in recommendation systems [31]. For the edge prediction task on real data, we also consider a supervised learning link-prediction routine (OLP) with topological predictors and the implementation of Ghasemian et al. [32] (see Appendix G 3 for details).

B. Performance for synthetic networks
We study various types of synthetic networks, generated by three different models to cover several network topologies. Two of them cover the extreme scenarios of networks generated, accounting only for community structure or only for reciprocity. For the former we use the standard stochastic block model (SBM) [9] and for the latter the reciprocity model of Holland and Leinhardt (HL) [10]. Our model, instead, is designed to tune the relative impact of community structure and reciprocity in determining edges, by varying the parameter η. Thus, we use the benchmark generative model described above to interpolate between these two extremes by tuning η: for small values we reproduce the results equivalent to the stochastic block model, whereas for higher values we replicate a structure similar to Holland and Leinhardt's model.
The generative process is described in detail in the Appendix A. As a remark, the exact joint likelihood of CRep is not determined in closed-form, however all the models used here for comparison adopt either its Poisson conditional distribution (our model with η > 0) or its Poisson marginal distribution (all the other models). Thus experiments here are aimed at highlighting differences in the various models' assumptions. By varying the network sparsity and the impact of communities and reciprocity, we illustrate types of structure that may exist in realworld data, and test each algorithm's robustness against them on various tasks including edge prediction and the ability to reproduce sample networks that replicate relevant network quantities.
Reproducing the topological properties An important property of a model is the ability to generate network samples that resemble what is observed in real data. We test this ability by considering topological properties like degree distribution, reciprocity, and hierarchical structure. We calculate their values on network samples, which are generated with the various generative models, by applying the inferred parameters from the given input data. Specifically, we consider networks generated synthetically as explained above, and for each individual network we infer the parameters by each model, and use them to generate five network samples. We compare topological properties of these samples with those observed on the ground truth networks used to infer the parameters.
In particular, we are interested in measuring reciprocity, as the networks generated by algorithms only based on community structure are not capable of reflecting the observed value of the reciprocity in the ground truth network, a shortcoming of these models which indeed limits their applications. The empirical evidence of this observation was part of the motivation to study this problem. In the experiments, we use the standard definition of reciprocity r, i.e., the ratio of the number of edges pointing in both directions to the total number of edges in the graph (we use the python implementation in networkx). As anticipated, in networks generated with the stochastic block model, r is often close to 0. Instead, a more interesting scenario is that of networks generated with the main purpose of replicating reciprocity, as in the HL model. This is an example of an exponential random graph model where reciprocity and sparsity are the two topological properties controlled in input. It is also one of the few cases where this type of model is analytical, see Appendix E. In this model, r is tuned by a parameter α so that the higher its value, the higher the reciprocity. Notice that, as usual in exponential random graphs models, latent variables such as communities are not considered. This model generates unweighted networks, hence r ≡ r w . Figure 2 shows that CRep significantly outperforms all the other generative models in reproducing r w , panel (a), and r, panel (b), as measured on the sampled networks.
The gap between the values of r and r w on the sampled networks is due to the mismatch between the binary adjacency matrices of the networks generated with the HL model (input data) and the weighted sampled ones generated with the various generative models, which use Poisson distributions. Similar results are obtained for the networks generated with our benchmark generative model. Also in this case, CRep captures reciprocity significantly better than the other models, consistently over a range of values of η as the input parameter. Moreover, in the case of fixed η, varying the sparsity and degree of overlapping communities lead to the same results. We leave details in the Appendix F 1.
At this point, we turn our attention to topological properties other than reciprocity, to investigate how these generative models perform in reproducing various relevant properties that might be of interest for a practitioner. Indeed, other possible mechanisms underlying network interactions are those that involve more than two individuals (which is the case for reciprocity), e.g., hierarchical structure, which requires the whole network for its computation.
As in our experiments we find that all models are able to retrieve the degree distribution with good accuracy, we mainly focus on replicating ranking of nodes, an application relevant when nodes have a score representing some intrinsic notion of relative strength or prestige. For this, we use SpringRank [33], an algorithm for inferring hierarchies in directed networks that assigns real-valued scores to nodes. We calculate the Gini index on these scores to provide a global measure for the whole network. Comparing the average over the five samples, we find that CRep and CRep0 are able to perfectly retrieve the Gini index of the original network, while the other models tend to overestimate it, see Appendix F 1. This is consistent over the various synthetic network topologies. Notice that this topological property is influenced neither by the value of η, nor the fraction of nodes with mixed-membership used to generate networks; however, it decreases as the average degree, and α increase.
Edge prediction We test the algorithms' ability in edge prediction tasks, in both cases of conditional and regular edge prediction. As we can see from Fig. 3, our model outperforms the others in conditional edge prediction, showing that it is able to efficiently exploit the additional information about the existence of the opposite edge. The performance gap between different approaches increases with η, as for high values of η, the reciprocity plays a bigger role in edge formation. In the opposite scenario of low η, the impact of reciprocity becomes negligible compared to community structure, and in this case we reproduce the same results as for the other algorithms. This is expected as our model infers small values of η in this case, thus in practice reducing to a conditional independent model as the others. Performance in terms of regular edge prediction is comparable to the other algorithms for small η, while it drops for intermediate values and then increases again as η grows. (b) Standard reciprocity r. Notice that r ≡ r w for the input data, but this is not true for the samples, as the generative models considered here generate weighted edges, i.e. the matrix A is in general not binary. Error bars are smaller than marker size. Unless otherwise stated, this will be the case in all of the figures.
These synthetic tests suggest that working with conditional probabilities results in more robust estimates of the probability that an edge exists if we have access to the edge in the opposite direction. Performance improvement is more significant when community structure is not the predominant mechanism in edge formation. We leave more details in the Appendix F 2.
To summarize results on synthetic networks, CRep is capable of suitably capturing the reciprocity values observed in a given network, while also retrieving hierarchical structures. Furthermore, CRep exploits the availability of extra information in performing edge prediction, by increased performance and robustness across various parameters' ranges.

FIG. 3: Edge prediction in benchmark networks.
Synthetic networks with N = 2100 nodes and K = 3 communities of equal-size unmixed group membership generated with the benchmark generative model proposed above by varying the reciprocity parameter η. The results are averages and standard deviations over three independent synthetic networks and over 5-fold of cross-validation test sets. The accuracy of edge prediction is measured with AUC and the baseline is the random value 0.5.

C. Performance for real networks
Above, we evaluated the ability of our model, CRep, to generate network samples that have reciprocity values as expected in input and tested its performance in edge prediction. In this section, we examine these abilities on real world datasets. We apply our method to datasets from a diverse set of fields, with sizes ranging up to N ∼ 10 4 nodes and up to E ∼ 10 5 links (see Table I and Appendix G 1 for details). Together, these examples cover various types of social relationships, communication interactions, transportation systems, and patterns of citations.
Reproducing the topological properties We apply the same procedure as before to infer the parameters Θ = (u, v, w, η) from data (this time, real networks) and then generate synthetic network samples based on them. Also in this case, CRep greatly outperforms the other models in reproducing r, consistently across datasets. We show as an example in Fig. 4 the results on the Erasmus dataset (Erasmus Mobility Network 2014 − 2018) [34], and we leave the others in the Appendix G 2.
Previously, we have discussed network-related quantities controlled by η, such as the expected fraction of edges purely due to communities (cr ratio ) or the quantity r w . Here we illustrate how the various real networks differ in the inferred values of η, which we denote asη. In particular, we show in Fig. 5 howη varies according to the reciprocity of these networks, unveiling a non-trivial pattern. While we see a general trend ofη increasing with r, there are interval ranges of r for whichη varies widely across networks, and vice-versa. For example, we see that for r ∈ [0.6, 0.8],η ranges in [0.1, 0.7]. This high variability suggests that r is the result of a complex combination of communities and reciprocity. We notice, for instance, that for high school friendship networks (HST and DT),η is low (i.e., in [0.1, 0.3]), showing that many reciprocated edges are explained by community structure. Instead, for online dating (POK) and communication networks (EU and DNC), we observe high values ofη, signaling a lower impact of communities, as reciprocity plays a bigger role. This reinforces the need to include in network models both mechanisms for explaining edge formation. Notice that these results are possible not only because our model accounts for reciprocity through an explicit parameter η, but also because it infers reciprocity values close to the observed ones, while the other methods fail at this, see Fig. 12.   Table I.
Edge prediction In the absence of ground truth, as in most real world networks, we test the ability in edge prediction by cross-validation, as done for synthetic networks. Table III shows the results in terms of AUC for the generative models CRep, MT, BPMF, as well as for OLP; the latter is a type of supervised learning technique which uses network topological information as features to predict the entries of A. CRep and OLP show the best results, with CRep having high performance for social networks. However, if we consider the conditional AUC, then CRep outperforms all the others in the majority of the datasets, as also observed in synthetic data. Finally, by averaging the AUC across the dataset, we find CRepnc is the best model. This confirms the ability of our model to efficiently exploit the additional information from the adjacency matrix to boost performance in terms of edge prediction.

IX. Case study: application of CRep to the Erasmus student exchange network
We illustrate our model on a real dataset to show various analysis that a practitioner can perform. We consider a network representation of the Erasmus student exchange program in 2018 [34], denoted as ERs18 in Table I. A node represents a higher education institution and an edge between nodes i and j denotes how many students were sent from i to spend a portion of their academic year abroad at institution j, as part of their study program towards a degree (Bachelor, Master, or PhD). This program is supported by the European Commission and involves N = 4389 institutions (mainly European), with a total of M = 90972 participating students in 2018. We recover community partitions from the network data using both CRepnc and MT, they have similar and high performance in edge prediction according to AUC (see Table III), and we fix K = 6 communities from crossvalidation. In Fig. 6, we notice that while both models find several groups that closely correlate with countries, CRepnc tends to put German institutions (left triangles) more in the same group (blue) and shifts few institutions in the red group, which seems made of mainly universities with strengths in engineering and technology (e.g., Universitat Politecnica de Catalunya, Politecnico di Milano and Institut Polytechnique de Grenoble). For instance, Università di Bologna, Federico II di Napoli and Padova have lower u i,red than what is predicted without accounting for reciprocity, instead Slovenská technická univerzita v Bratislave, Kauno Technologijos Universitetas and Universidad de Oviedo increase their membership in this group.
In addition, CRepnc places more institutions with higher membership in the green group, see Fig. 6 (g) (hard membership). While there is no apparent common attribute between these (e.g., country), we find that many nodes with high "green" entry of u i tend to reciprocate more edges. Specifically, they have a high fraction of out-neighbors such that λ 0 ij is much smaller than λ 0 ji . That is, the edges A ij such that A ji also exists, have a lower impact in determining the value of u i in the algorithm. In fact u ik ∝ j,q A ij ρ ij +η Aji , see Eq. (7). Hence, if the denominator is high because of A ji , the weight of the edge A ij decreases. Nodes with many such A ij tend to have lower entries u ik and thus lower λ 0 ij . This is a qualitative explanation for having different membership, however the situation is more complicated than this, as one needs to account for the effects on the whole network. In fact, also v jq changes between the two algorithms, for a similar reason, thus also contributing to a different u ik .
The primary benefit of CRep, however, lies not in its ability to recover the communities but in what it reveals about the reciprocity patterns in the network. Home and receiving institutions must sign an inter-institutional agreement to allow for student exchanges between them. While institutions may sign them because of clear affinities between their educational training offerings (e.g., both universities are strong in natural science), they might also do so because of some mechanisms involving reciprocity, as hosting students costs resources. Moreover, reciprocity could be further increased by previous knowledge or collaborations between individual faculties, thus institutional reciprocity may be also driven by faculty reciprocity. In addition to the communities themselves, our model also returns η, which can reveal features of the data related to such reciprocity effects not seen with standard generative models, such has cr ratio or E [A ij |A ji , Θ]. We find a maximum likelihood value of η = 0.4, signaling a significant reciprocity effect. In fact, according to Eq. (14), on average 40% of the edges are influenced by reciprocity.
While η gives a global picture of the whole network, our models still allows to distinguish the impact of reciprocity on individual edges. For instance, if an institution i accepts many students from j, then j might be more willing to accept students from i, even though i's features might not match j's preferences. If we distinguish the u i as the set of preferences of i and v j as the set of attributes of j, then our model will naturally convey this through high λ 0 ij and low λ 0 ji for such a case. CRep is able to capture these situations quantitatively, by means of the quantities cr ij := λ 0 ij /m ij (a cr ratio per edge) with values in [0, 1] which measures the relative contribution of communities alone to determine edges between i and j. Focusing on a single institution i, one can analyze the difference d ij := cr ij − cr ji ∈ [−1, 1] for all j such that both A ij , A ji > 0 and find different reciprocity patterns, as we show in Fig. 7. Here we plot three extreme cases where i has most of the d ij being less, equal or greater than 0. The Universidad Pablo de Olavide in Sevilla, panel (a), has mostly d ij < 0 (plotted in red), meaning that reciprocity has a strong effect in determining its out-going edges to universities that instead send students to Sevilla mostly out of community preference. The opposite case is that of Technische Universität München, panel (b), which has most of the d ij > 0 (plotted in blue), signaling that it tends to select its out-going edges more out of preference than their counterparts, who tend to reciprocate instead. Università degli Studi di Firenze, panel (c), is an example of an institution with several d ij close to 0 (plotted in white), meaning that most of its reciprocated edges are due to community affinities. In other words, Firenze selects out-going j based on preference and those who select Firenze do the same, so the impact of reciprocity is low. Apart from these three extremes, many universities display a range of such behaviors; we give an example of Universidad Carlos III de Madrid, panel (d), which has a balanced fraction of reciprocated edges covering these three cases (there are about 1/3 of blue, red, and white edges in the corresponding figure). Notice that the value of d ij yields an incomplete picture of the situation, since it does not distinguish between cases where the quantities cr ij , cr ji have different magnitudes while keeping their difference constant.

X. Conclusion
CRep is a mathematically principled generative model for capturing both community and reciprocity patterns in directed networks. It relies on relaxing strict conditional independence assumptions on edges that limit the applicability of standard methods on real problems where reciprocity plays an important role. Its algorithmic implementation is efficient and scalable to large system sizes. The corresponding generative model allows for the creation of synthetic networks with the desired interplay FIG. 6: Erasmus 2018 community structure. For visualization clarity, we show the subnetwork made of the 10% biggest institutions and the 3000 edges with highest weights (inference was performed on the whole network). Panels (a)-(f) show groups K = 1, 2, ..., 6 (mixed membership) by CRepnc, and panels (h)-(m) show the same groups by MT. Panel (g) illustrates the groups by CRepnc in the case of hard membership, while the groups by MT are represented in panel (n). Node color intensity increases with u ik , so that darker nodes have stronger membership u in that group, each color is a group (mixed membership) and nodes with light blue border are nodes that change the most the membership in the two algorithms; for each group k, we only show nodes that have u ik > 0.1. Node and edge size are proportional to the size of an institution measured by the total number of outgoing and incoming students. Node shapes denote country. between community and reciprocity in determining the edges, while allowing the tuning of network sparsity.
In addition to providing all the analysis tools typical of standard generative models with communities, our model makes it possible to answer questions about reciprocity in networks that were not previously possible; for instance, performing probabilistic conditional edge prediction and estimating the relative contribution of community and reciprocity in determining edges. We show how real networks display a wide range of the reciprocity parameter, signaling the variety of possible patterns for this property. In the context of the Erasmus student exchange network, our model allowed us to distinguish universities based on their pattern of reciprocated edges.
More generally, our model shows how we can relax strict conditional independence assumptions on edges and showcases possible consequences in doing this. This presents an opportunity for researchers to rethink the fundamental assumptions behind generative models, and present models that may open doors to new theories and questions. We make one step in this direction, as our model connects two popular problems that are mainly treated independently: the inference of communities in networks and generating directed networks where reciprocity plays a relevant role. We used this connection to obtain networks with community structure and values of reciprocity consistent with those observed in real data.
Both the assumption and the model we have presented are only the first step in a broader line of work that investigates how certain topological properties are reflected in networks with latent community structure as dominant mechanism in edge formation. There are a number of directions in which this work could be extended. We have considered here a simple way to account for reciprocity and break conditional independence, by considering a unique parameter for the whole network. Our model could be extended to account for node-dependent parameters, where reciprocity varies between individuals. In addition, possible extensions may incorporate extra information such as degree, attribute or signals on nodes [30,[35][36][37][38], edges of different types as in multilayer networks [16] and dynamics in time [39][40][41][42][43][44]. Reciprocity is one of the many effects that could play a role in determining how nodes interact in a network. One could go further Node size is proportional to university size, the shape denotes country, the colors are the highest entry of u i (for the four reference nodes -white node border) and v j (for all its neighbors). Edge size is proportional to its weight; edge colors vary continuously from red to blue, based on the value of d ij = cr ij − cr ji : high intensity red, white and high intensity blue mean close to -1, 0 and 1 respectively. than this by considering incorporating quantities that account for triples of individuals, for instance clustering coefficient, transitivity or global centrality measures [45]. These properties cannot be captured by standard SBMlike models [46]. In this respect, a recent work of Peixoto [47] shares some similarities with ours considering triadic closure instead of reciprocity, making an effort towards extending the stochastic block model framework to incorporate more elaborate topological structure that is not captured otherwise. This is something that exponential random graphs or stochastic actor oriented models are capable of [14,[48][49][50][51], without including latent community structure but rather fitting network statistics. In probabilistic generative models, this would require further breaking conditional dependencies between edges, potentially increasing the model complexity to encompass more complicated situations. With our work, we made the first step in this direction.
While there is no unique generative model that captures all the possible network properties well, our work illustrates how to target reciprocity. As our original motivation to study this problem came from the realization that standard generative models fail to generate synthetic networks with meaningful values for this property, our work illustrates a way in which latent variable frameworks can be applied more realistically, and provides an example of how network scientists can better align fundamental theories with realistic applications. We provide an open source implementation of the code online at https://github.com/mcontisc/CRep.
The synthetic networks used in the analysis are of three types and represent different scenarios: networks with community structure only, with reciprocity only and networks with both communities and reciprocity. In order to obtain networks with only a community structure we use a stochastic block model with different values of average degree k . We generate networks with K = 3 communities of equal-size unmixed group membership, N = 2100 nodes and an assortative structure (w has higher diagonal entries) with main probabilities p 1 = c K/N and entries outside the main diagonal equal to p 2 = 0.1 p1, so that the average degree is k = c+(K −1) c/10, where c is the average degree within the same community. We generate three independent samples for each value of c ∈ [2,20], that corresponds to k ∈ [2.4, 24]. On the other hand, we generate networks influenced by reciprocity only through an implementation of the reciprocity model proposed by Holland and Leinhardt (see Appendix E for details). The input parameter α can be tuned to obtain different values of network reciprocity and we generate three independent samples for each value of α ∈ [0, 10]. We consider N = 1000 nodes and a probability to generate one of the directed-edges equal to p = 0.002.
In order to work with synthetic networks having an intrinsic community structure and a given reciprocity value, we use the benchmark generative model proposed in this paper. We generate networks with N = 2100 nodes and K = 3 communities by varying three different input parameters: the average degree k ∈ [2,20], the reciprocity coefficient η ∈ [0, 1) and the fraction of nodes with mixed membership over ∈ [0, 1]. While varying one of the parameter, the others are fixed to k = 20, η = 0.5 and the degree of overlapping communities over = 0. In detail, networks are generated in two steps. First, membership vectors u and v are generated following an equalsize unmixed group membership and a Dirichlet distribution with parameter α = 0.1 for the entries with mixed membership; and the affinity matrix w is generated using an assortative block structures with main probabilities p 1 = K/N and secondary probabilities p 2 = 0.1 p 1 . Thus the latent variables Θ = (u, v, w, η) are fixed. Second, edges are drawn according to the generative model described in the main text. Specifically, for each pair of nodes (i, j), i) extract A ij from a Poisson of mean as in Eq. (10); ii) extract A ji from a Poisson of mean as in Eq.
(3). This procedure results in a directed network with the desired reciprocity and sparsity. We generate three independent networks for each value of the three different input parameters.

B. Edge prediction and cross-validation
We perform edge prediction using 5-fold crossvalidation. In each realization, we divide the dataset, i.e., the entries A ij of the adjacency matrix, into five equal groups selected at random. We use four of these groups as a training set, to infer the parameters Θ. We then use the fifth group as a test set, evaluating the score for each A ij in this set, and calculate the AUC value. By varying which group we use as the test set, we get 5 trials per realization. The final AUC is the average over these. To compute the regular AUC we use as score the expected value E P (Aij |Θ) [A ij ] = m ij as in Eq. (10); for the conditional AUC (AUC−cond), we use as score E P (Aij |Aji,Θ) [A ij ] = λ 0 ij + η A ji , i.e., the expected value over the conditional distribution. Notice that the latter can only be computed for CRep, as for the others m ij ≡ λ 0 ij , and thus the two AUC values coincide. The AUC is specified for binary entries, thus the edge weight is not accounted in the evaluation. However, our goal here is to assess edge existence, hence AUC is a suitable metric for this. If a practitioner aims at assessing the quality of the inferred weights as well, then one should specify different metrics for this.

C. Inference: numerical implementation
All the generative models require inferring K, the number of communities. We select this by cross-validation. Specifically, we run several held-out trials as explained above by varying K and select the value of K that gives the highest (regular) average AUC on the test sets. We then extract the parameters of each method using their best K. For MT, BPMF and CRep0, we extract the parameters u, v, w; in addition, for CRep and CRepnc, we extract η. All these algorithms converge to a local optima, as the likelihood landscape is not convex. Hence, we run the algorithm 10 times for different random initializations of the parameters and select the realization that has higher likelihood value.

D. Detailed derivations
We derive in detail the equations for inferring the parameters. We first apply a variational approach to make the problem tractable, and then use an expectationmaximization algorithm to derive the equations of the updates.

Variational approach
We aim at maximizing the log-pseudo-likelihood in Eq. (5). The first step is to facilitate the maximization process of the logarithmic term. We consider a probability distribution ρ ij over the two competing terms: this is our estimate of the probability that the edges exist due to the contribution of either the community membership or the reciprocity term. Applying Jensen's inequality logx ≥ log x: (D1) Moreover, this holds with equality when: Thus maximizing L ps (Θ) is equivalent to maximizing: We apply once more the variational approach to make the sum inside the logarithm tractable. Similarly as before, we introduce a probability distribution φ ijkq such that: The equality holds when: Thus maximizing L ps (Θ, ρ) is equivalent to maximizing: with respect to Θ, ρ, φ.

Expectation-Maximization updates
Equations for the updates of each of the parameters can be obtained by taking the derivative of Eq. (D5) with respect to a given parameter and setting it to zero. For instance, the update equation for η is obtained by considering the partial derivative: Setting this to zero and defining M = i,j A ij , we obtain: Similarly, for the community affinity matrix we get: Here we show how to enforce constraints like k u ik = 1, which is an arbitrary choice that can be easily incorporated into our model. To this end, it is convenient to rewrite the log-pseudo-likelihood as follow, Then, following the approach in [52], to simplify the maximization of the log-pseudo-likelihood, we substitute w kq from Eq. (D8) into Eq. (D10): The second term in the above equation does not depend explicitly on u ik and v jq . In order to apply the constraint on the maximization, we add Lagrange multipliers γ u i , γ v i : The update equation for u ik is obtained by considering the partial derivative, and setting it to zero, which yields: By applying the normalization constraint on the u ik , i.e., k u ik = 1, and noticing that ρ ij +η Aji , we can find an expression for γ u i : Similarly, we have the following update equation for v: where Deriving the expected value of the marginal distribution Solving for m ij yields: which implies: With similar calculations as before we obtain: To fully determine this expression we need to specify the second moment E A 2 ji . For binary variables, we could assume E A 2 ji = E [A ji ] = m ji , as this is the case for Bernoulli distributions. With this assumption, we obtain E [A ij A ji ] = λ 0 ij + η m ji . Alternatively, we can assume E A 2 ji = m ji +m 2 ji as is the case for the Poisson distribution, and thus obtain E [A ij A ji ] = λ 0 ij + η m ji + η m 2 ji . Finally we have: where in the first row we use the first order Taylor expansion as an approximation. With this assumption, we obtain that the parameter η is a lower bound for the expected value of r w . An equivalent expression can be derived for models that assume conditional independence, e.g., our model with η = 0. In this case we get: which yields:

E. Holland and Leinhardt reciprocity model
The model assumes an unweighted and directed network, i.e., asymmetric adjacency matrix with binary values A ij ∈ {0, 1}, and the following joint probability: where Z(θ, α) = 1 + 2e −θ + e −2θ+α is the normalization term. The parameter α controls the level of reciprocity, it couples the two entries A ij and A ji thus making the model not factorized; edges between different pairs (i, j) are conditionally independent given the parameters. This is one of the few analytically tractable exponential random graph models. Due to this property, we can extract analytical marginal and conditional distributions for a pair of nodes (i, j): These expressions can be used to sample networks with the joint distribution given in Eq. (E2). Tuning the value of the parameter α, one generates networks with different values of reciprocity.
F. Performance in synthetic networks 1.
Reproducing the topological properties Here we show in more details the ability of the models to reproduce network samples that replicate relevant network quantities. Figure 8 shows r and r w as defined in Eq. (15), computed in the sampled networks of synthetic data generated with a stochastic block model and our benchmark generative model. As expected, the reciprocity in networks generated with the stochastic block model is always close to zero. Instead, the networks generated with our benchmark generative model present different values of reciprocity, and CRep captures these values significantly better than the other models, consistently across various magnitudes of input η. Even in the case of fixed η, by changing sparsity, we observe the same pattern. By varying the degree of overlapping communities we obtain the same results as changing the average degree (we do not report them here). Figure 9 shows the Gini index computed on nodes scores obtained with the SpringRank algorithm. The Gini index provides a global measure for the whole network, the higher its value, the more hierarchical the network is. We compare the average over the five samples, and we find that CRep and CRep0 have reasonable accuracy in retrieving the Gini index of the original network, while the other models tend to overestimate it. This is consistent over the various synthetic network topologies, i.e., network generated with the stochastic block model, panel (a), the HL model, panel (b), and our benchmark generative model, panel (c). Furthermore, we notice that this topological property decreases as the average degree within the same community, c, and α increase, while it is not influenced by the value of η. We omit the results for the networks generated with our benchmark generative model by varying the sparsity and the fraction of nodes with mixed-membership because we obtain similar results to the stochastic block networks and the benchmark data by varying η, respectively.

Edge prediction in synthetic networks
Here we show the results in terms of edge prediction on synthetic data generated with our benchmark generative model by varying the average degree k and the fraction of nodes with mixed-membership, which we denote over. We use both conditional and regular edge prediction and Fig. 10 highlights the robustness of CRep and CRepnc in terms of conditional edge predictions, as their performance are significantly higher than that of the other algorithms and do not decrease with increasing overlapping communities and sparsity. Indeed, the results are robust, as we vary the fraction of nodes with overlapping community membership and the average degree, while fixing η = 0.5. Notice also the stability of CRep and CRepnc in terms of regular edge prediction and how they outperform the other models in critical ranges, e.g., small k and high over.
Moreover, we find more stable results also in terms of regular edge prediction, where CRep and CRepnc have constant values across the different input parameters, outperforming other methods in critical ranges, e.g., small average degree or high overlap between communities. The results of our experiments suggest that working with conditional probabilities results in more robust estimates of the probability that an edge exists if we have access to the edge in the opposite direction. Performance improvement is more significant when community structure is not the predominant mechanism in edge formation.

Community detection in synthetic data
For sake of completeness, here we show the performance of the models on recovering communities. We consider as performance measure the F1-score (F1) and cosine similarity (CS), the former one is valid for hard membership while the latter captures mixed-membership, we calculate for both the average over the nodes. When measuring the F1-score we consider the entries of maximum value of the membership vectors. Both measures are between 0 and 1 and a value of 1 means perfect reconstruction. Figure 11 shows the accuracy in networks generated with the benchmark generative model by varying the reciprocity parameter η and for synthetic data created with a stochastic block model by varying the average degree within the same community c. For comparison in these last networks, we consider also the Leiden algorithm [53], a non-generative method. Even if community detection is not the main focus of our model, we notice the ability of CRep in retrieving communities in networks without reciprocity, while its performance decreases as reciprocity increases. This is expected as the community impact in determining the likelihood of an edge decreases as η increases. Notice that the benchmark data have been generated with fixed k = 20, thus models without reciprocity are capable of fully recover- Synthetic networks with N = 2100 nodes and K = 3 communities of equal-size unmixed group membership generated with the benchmark generative model proposed above by varying (a) the average degree k and (b) the fraction of nodes with mixed-membership over. The results are averages and standard deviations over three independent synthetic networks and over 5-fold cross-validation test sets. The accuracy of edge prediction is measured with AUC and the baseline is the random value 0.5.
ing the community even in the case where reciprocity is there, provided that the average degree is large enough. These synthetic tests suggest, on one side, the robustness of community detection-only methods in recovering communities even in the presence of reciprocity; on the other side the good performance of CRep in recovering communities when reciprocity has intermediate or low level. This is somehow expected, as this model gives increasingly less weight to the communities as reciprocity increases, thus it is not optimized to recover the communities when these are not fully determining edge formation. We apply our approach to different types of networks, such as social, infrastructure, online communication, and citation networks. Table I provides a brief overview of the datasets studied in this work, as well as their ab-breviations. All datasets, have been pre-processed as follows: i) self-loops are removed; ii) only nodes that have at least one out-going and one in-coming edge are kept; iii) we used only the giant connected components. Some datasets require additional specific preprocessing. Specifically, the citation networks (here: CIT05, SCC2016, ACMv9 ) require extracting a network author-author from a network of paper-citation, so that an edge means that an author cites another author. Furthermore, we split dynamic networks into separate individual networks where we kept only interactions happening within a certain time window. This applies to Dutch (DT2, DT6), High school friendships (HST11, HST12, HST2), online dating (POK0, POK6, POK12), and Erasmus (ERs14, ERs15, ERs16, ERs17, ERs18).

Reproducing the topological properties
Here we show the ability of the models to reproduce network samples that replicate relevant network quantities. For each real network we infer the parameters by each model, and use them to generate five synthetic network samples. Figure 12 shows the reciprocity r. For each model, it outputs the averages and the standard deviations over the five samples and the dashed red lines indicate the r value of the input datasets. We notice the heterogeneity of the analysed networks and how CRep adapts to all different situations, while the other models underestimate the true value most of the times. Figure 13 shows the Gini index computed on nodes scores obtained with the SpringRank algorithm. The results vary widely depending on the datasets, and we cannot draw general conclusions. In this scenario, we have also studied the reproducibility of the clustering coefficient, i.e., the tendency of nodes to form edges within the same neighborhood, however, we obtain poor results in line with the SBM approach, as predicted in [46]. Moreover, these are topological properties that involve more complex interactions than pairwise, as in the case of reciprocity (clustering involves triangles and SpringRank score is a global measure). This suggests that, in order to have better performance, one would need to develop more complex models, for instance extending the ideas behind CRep to capture triadic interactions, possibly guided by domain-knowledge about how triadic interactions and reciprocity are related [45]. We leave this for future work, noting that while exponential random graph models can do this, they do not include latent community structure (analogously as for reciprocity).

Link prediction features
Here we present the supervised learning-link prediction routine (OLP) used for comparison in the edge prediction task on real data. In the link prediction task, scores are assigned to all possible pairs of nodes in the graph based on a set of criteria. Then, the pairs of nodes are sorted according to their scores in an ascending order and the most-likely links are the pairs with scores above a threshold value.
Two categories of features are used to determine the criteria of link classification: (i) global features, defined based on the features of the entire network, such as the number of nodes, number of edges, average degree of nodes, and the average clustering coefficient, and (ii) local features, which include the descriptive features of a single node or a pair of nodes.
In this work, we apply the extended definition of features for a directed network of Ghasemian et al. [32]. We also examine the effect of belonging to the same community on the local pairwise features, i.e., pairwise attributes contribute in the link prediction only if the two nodes belong to the same community. However, we did not find significant changes and at the price of higher computational cost, hence, we exclude this factor from the study and omit the results. Considering Γ (x) out/in as the set of out/in-neighbors of node x, and d(x, y) as the distance between nodes x and y, some of the wellknown features deployed for link prediction are presented in Table II. FIG. 12: Reciprocity in real networks. Empirical averages and standard deviations of reciprocity r over 5 samples of each real network (see Table I for details). The red dashed lines indicate the r on the input networks.  Table I for details). The red dashed lines indicate the values on the input networks.

Resource Allocation index
defined for a pair of nodes: x, y: z∈{Γ (x) out/in ∩Γ (y) out/in } 1 |Γ (z)| Betweenness centrality a measure of node centrality based on the shortest paths Closeness centrality defined for a pair of nodes: x, y:

Shortest Paths
shortest path between nodes: x, y Katz centralities a measure of centrality in a network PageRank centralities a measure of the importance of a node as an adjustment of Katz centrality Eigenvector centralities an adjustment of Katz centrality of a node in regards to the importance of its neighbors Clustering coefficient for node x number of triangles connected to node x number of triples centered around node x Preferential attachment the tendency of nodes to connect to the nodes with higher degree Common community 1 if the pair of nodes belong to the same community, otherwise zero  Table I for details). Results are averages and standard deviations over 5-fold cross-validation test sets.
In grey box we show the best performance over all methods, while in boldface the best results in terms of regular AUC. The last row reports the average and standard deviation of each method over datasets.