Identifying maximal sets of significantly interacting nodes in higher-order networks

We introduce a method for the detection of Statistically Validated Simplices in higher-order networks. Statistically validated simplices represent the maximal sets of nodes of any size that consistently interact collectively and do not include co-interacting nodes that appears only occasionally. Using properly designed higher-order benchmarks, we show that our approach is highly effective in systems where the maximal sets are likely to be diluted into interactions of larger sizes that include occasional participants. By applying our method to two real world datasets, we also show how it allows to detect simplices whose nodes are characterized by significant levels of similarity, providing new insights on the generative processes of real world higher-order networks.

Over the last decades we have witnessed an explosion in the number and size of data describing the different domains of the world we live in. In particular networks, collection of nodes and edges linking pairs of them, have emerged as a very valuable tool to describe and model the topology of all such real-world systems where interactions among agents play a key role, from trade economics to complex biological processes [1]. Nevertheless, mounting evidence is now suggesting that networks can only provide a limited description of many complex systems, as they are inherently limited to model pairwise interactions [2,3]. Higher-order interactions have been observed in a variety of contexts, from human interactions in face-to-face systems [4] to collaboration networks [5], ecological systems [6] and the human brain [7,8]. Interestingly, such non-dyadic interactions were shown to give rise to new collective phenomena in dynamics such as contagion [9,10] and cooperative process [11], as well as diffusion [12,13] and synchronization [14][15][16].
To properly encode such higher-order systems [17], more complicated mathematical frameworks are needed, such as simplicial complexes and hypergraphs. In hypergraphs, generalised hyperedges encode interactions among an arbitrary number of elementary units. In simplicial complexes higher-order interactions are described by simplices that satisfy a particular inclusion property, such that also all subfaces of a given simplex must exist. Practically, this means that if n interacting nodes are described by a simplex, also all interactions among all possible subsets of nodes are considered to exist. This additional constrain comes at the benefit of enabling the use of computational techniques from algebraic topology, typically known as topological data analysis, to efficiently describe the topological features of such systems [18,19]. Given the growing interest in higher-order interactions, over the years several network metrics have been generalised to account for the presence of non-dyadic interactions. These include local scores such as centrality measures [20][21][22], higher-order motif analysis [23,24] and simplicial closure [25], the proposal of methods to extract higher-order community structure [26][27][28][29], or unveiling higher-order dependencies in multivariate time-series [30].
The rising interest towards higher-order complex systems calls for the design of suitable filtering techniques. In fact, in light of the explosion of data on systems belonging to different domains, filtering methods have attracted attention as tools that allow to focus on relevant connections that are highly representative of a system, while discarding all the redundant and/or noisy information. This approach is useful to detect relevant trends in the data, otherwise hidden by noise or by the intrinsic heterogeneity of investigated complex systems.
Filtering approaches for complex systems can be divided into two main families. The first family of filtering approaches relies on imposing topological constraints that force to reduce the number of connections in the filtered system. The first example of such approaches is the Minimum Spanning Tree, which reduce a system of n connected nodes to a tree with n − 1 links that are placed starting from those with the highest weight [31]. This approach has been extended by proposing new topological constraints [32,33] and the robustness of the filtering procedures has been assessed through statistical procedures [34,35].
In the second family we find methods that drop the connections that are compatible with an appropriate null model. Null models are used to produce ensembles of replicas of the system, where some key feature are preserved while the generative model that produce links and/or hyperlinks is randomized. An example is the configuration model, in which the degree of each node is exactly preserved but the edges are placed randomly. Different null models have been proposed in the last decade, for different classes of systems. There are null models for weighted networks [36][37][38], bipartite networks [39][40][41][42][43][44], temporal networks [45] and correlation matrices [46,47]. Recently, the first null models for higher-order systems have been proposed [45,48].
Statistically validated filtering of complex systems has been applied to several domains, including actor/movies networks [39], trading decisions of investors [49][50][51][52], financial transactions between major market agents [53], mobile phone calls of a large set of users [40,54], financial credit transactions occurring in an Interbank market [41], intraday lead-lag relationships of returns of financial assets traded in major financial markets [55], the Japanese credit market [56], features of communities detected in complex networks [57],international trade networks [58], social networks of news consumption [59], and rating networks of e-commerce platforms [60].
In this work, we propose a filtering methodology for higher-order complex systems that is able of extracting Statistically Validated Simplices from collectively interacting groups of nodes. Specifically, we develop a method using a null model designed for repeated higher-order interactions, that is able to identify the maximal sets of nodes that consistently interact collectively and discard the nodes that co-interact with them only occasionally. We label these minimal sets of nodes as Statistically Validated Simplices, as the filtered significant units share the inclusion property that defines a simplex. Since our null models expands a class of existing models [39,48], it inherits from them the capability of preserving the heterogeneity of the system. The paper is organized as follows: (i) in section I we introduce and motivate our method and discuss a possible approximation of the proposed p-value; (ii) in section II we introduce a class of synthetic benchmarks that we use to test the accuracy and precision of our method against available alternative in literature. We also apply our method to empirical data and show that it is able of extracting relevant information on the investigated systems; (iii) finally in section III we draw our conclusions.

I. FILTERING METHOD
We aim at obtaining a tool that is capable of identifying groups of agents that are consistently interacting among them. Specifically, we want our groups to contain agents that are consistently appearing together in collective interactions, without including those agents that might appear in some of these interactions but whose presence is occasional. For example, consider a set of n products that are consistently bought together in the same basket because of functional relatedness. Given the heterogeneity of shopping habits, it is likely that these baskets will feature occasional other products. We want our tool to identify the n products that constitute a cohesive unit and discard all those other ones that are copurchased only occasionally. In the following we adopt the formalism of higher-order networks [2], that naturally provides the concepts required to build our tool. Within this formalism, agents are mapped into nodes and higherorder interactions between them are mapped into hyperlinks. We also recall the concept of simplicial complex as a collection of simplices, where a simplex of order n is a set of n + 1 interacting nodes, with the condition that all interactions between any subset of the nodes also exist. In what follows, to avoid confusion, we will always characterize hyperlinks and/or simplices through their size, which is the actual number of nodes that are interacting. A hyperlink/simplex of size n corresponds to an interaction between n nodes. Thus in this framework our research question can be rephrased into the problem of transforming a collection of hyperlinks that encode recorded higher-order interactions between a set of real world agents into a simplicial complex where each simplex represent a set of statistically validated interacting nodes. This transformation comes with significant implications. Indeed, a simplicial complex is a mathematical tool that allows to analytically generalize to higherorder systems many results that hold for pairwise networks [9,61]. However, not all real world collections of hyperlinks constitute a simplicial complex, as a simplicial complex K is valid only if, for every simplex of size n that constitutes it, all its sub-simplices of any size n < n also belong to K. A collection of real world hyperlinks in general does not hold this property, as there is no guarantee that for each set of n nodes interacting in a hyperlink all subsets of n < n nodes will also be interacting. On the contrary, we want our statistically validated groups of co-interacting nodes to be simplices because we consider all subgroups of n nodes that are significantly interacting together to be also significantly interacting. To keep our example, if n products are consistently bought together, any subgroup of these products will also be significant, as they will appear on the same baskets. In other words, we want to identify the statistically significant higher-order behavior of agents -i.e. the core of n products that are bought together consistently -, going beyond the detection of higher-order observables -i.e. all the specific combinations of core products and occasional ones that are sometimes bought together without constituting a significant basket (Fig. 1). In what follows, we will call the filtered set of significant sets of interacting nodes Statistically Validated Simplices (SVS). Given n nodes that are active in a system, we need a statistical criterion to assess whether the number of times in which they interact together is statistically significant. The simplest way to obtain this result would be to set a threshold on the number of times required to consider a group significant. However, this approach is subopti-mal because it is incapable of taking into account the heterogeneity in activity among agents. In fact, as the individual activity of agents in complex systems is usually strongly heterogeneous, a single threshold -e.g. one extracted from the distribution of activity of all nodeswould be too small for hyperlinks that contains the most active nodes and too large for hyperlinks that contains the least active nodes. Also, it would require to design an arbitrary set of rules to deal with hyperlinks whose participating nodes have heterogeneous levels of activity. For this reason, the validation of networks and hypergraphs is usually performed by means of a null hypothesis that is able to naturally incorporate the heterogeneous activity of nodes. Here we introduce a modified version of the null hypothesis that was designed for validating hyperlinks for Statistically Validated Hypergraphs (SVH) [48]. In fact, the SVH approach is limited to validate higher-order observables of a specific size n, i.e. observed hyperlinks of n nodes. The statistical significance of a hyperlink is assessed by considering its weight X -i.e. the number of times in which it is observed -, the number of hyperlinks of size n to which each node participates and the global number of hyperlinks of size n. From these parameters one can extract a p-value that represents the probability of observing the same weight X or more under a null hypothesis in which each node randomly selects its n − 1 partners. Here, we want to modify this approach by considering and testing all possible combinations of nodes that co-interact -no matter with how many other nodes -at least once. To obtain this, we modify the null hypothesis used in the SVH in the following way.
For simplicity, we start from the case of 3 nodes i, j, k that co-interact N ijk times in hyperlinks of size n >= 3. The three nodes appear respectively in N i , N j , N k hyperlinks. As we are not limiting anymore the possible hyperlinks to those whose size is exactly 3, we now set as N the total number of hyperlinks observed in the system. This choice for the value of N consider all hyperlinks equally likely to contain any number of nodes. Alternatively, N can be set as the number of hyperlinks with size larger or equal n -in our case n = 3. The hyperlinks that are counted in N ijk all contains the nodes i, j, k but, differently from the case of SVH, are not limited to them. For example, a hyperlink that contains the nodes i, j, k, l also contributes to N ijk . Under the null hypothesis that each node selects randomly the hyperlinks to which it participates -and thus its n − 1 counterparts in any hyperlinkthe probability of observing i, j, k interacting N ijk times is where H(N AB |N, N A , N B ) is the hypergeometric distribution that computes the probability of having an intersection of size N AB between two sets A and B of size N A and N B given N total elements. The probability p(N ijk ) in Eq. 1 represents the probability of having a random intersection of size N ijk between the three sets of hyperlinks of nodes i, j, k out of N total hyperlinks [62], and is obtained through the convolution of two instances of the hypergeometric distribution. Indeed, we start from the probability of having an intersection of size X between nodes i and j and multiply it with the probability of having an intersection of size N ijk between the hyperlinks of node k and the intersection set of size X between the hyperlinks of node i and those of node j. This product is then summed over all possible values of X, i.e., all possible intersections between i and j which are compatible with the observed number of interactions between all the three nodes, which lie in the interval Eq. 1 is symmetrical around permutations of i, j, k as the specific order with which they are considered does not affect the value of the probability. Starting from Eq. 1, we then compute a p-value for the symplex that contains i,j and k through the survival function, The p-value provides the probability of observing under our null hypothesis N ijk or more hyperlinks that contain -but are not limited to -the nodes i, j, k. For a generic simplex of n nodes, Eq. 1 becomes ... × In order to assess the statistical significance of a simplex of size n, its p-value is tested against a threshold of statistical significance α, after adding a multiple hypothesis test correction which is needed because of the high number of tests -one per each simplex. In all the results presented in this paper we use α = 0.01. However, since we consider all smaller combinations of nodes constituting a significant simplex to be themselves significant simplices, we proceed to p-value testing in a sorted order. In fact, we start from testing the largest simplex and we then proceed towards the smallest. If a simplex of size n passes our statistical test and is thus selected as informative because it rejects the null hypothesis, we do not test any of its smaller subsets. Moreover, as the number of possible simplices of size n from N nodes nodes explodes as N nodes n , we only test simplices whose n nodes appear in at least one observed hyperlinks. This means that we consider candidates for being significant simplices only sets of nodes that interact at least once. Thus, for each size n, we proceed from largest to smallest (where the largest n corresponds to the size of the largest observed hyperlink), and we test all candidates size by size after correcting the significance threshold by means of the control for the False Discovery Rate [63]. The total number of test considered at each size n is N t = N nodes n which is the number of all possible combinations with no repetitions of size n of the N nodes elements that are active in the system. Thus, when applying the control for the False Discovery Rate method we start from a Bonferroni threshold computed as α B = α/N t .

A. Approximated p-values
Extracting the probability of Eq. 3 requires the convolution of n − 1 instances of hypergeometric distributions, each of which is built through the combination of several binomial coefficients. For this reason, this task becomes computationally intense for larger n. On top of this, our approach may require the computation of several p-values over the same nodes, as any set of n nodes which does not represent a significant simplex is split into n candidate simplices of size n − 1. If not significant, each of these complexes is further split into candidate simplices of size n − 2, and so on, quickly causing the number of computationally intense p-value that need to be computed to explode. To overcome this computational bottleneck, we have tested an approximated version of the p-value. Eq. 3 represents the probability that n nodes appear in N 12...n hyperlinks after each of them has randomly sampled its hyperlinks from an urn. As the random sample performed by each node is independent from the others, we can approximate this estimation using the binomial distribution, turning Eq. 3 into Ni N . In other words, we are estimating the probability of observing N 12...n hyperlinks (out of the N total hyperlinks observed) that involve the n nodes, where each hyperlink has a probability of being observed which is the product of the normalized number of hyperlinks Ni N in which each node appears. Our approximated p-value recalls the one proposed in Ref. [45], as both approaches use the binomial distribution. However, here we model differently the probability of observing a collective interaction. In Ref. [45], given the specific nature of the investigated data that only tracks temporal pairwise interactions, the p in the binomial distribution is obtained by multiplying the probabilities of observing all pairwise links that constitute a clique of the same size of the tested hyperlink. Here, we obtain it by multiplying the normalized number of hyperlinks of each node involved in the interaction. The approximation of Eq. 4 lies in the fact that we do not take into account anymore the n-dimensional array of individual activities of each node in the hyperlink {N i }, but we consider the normalized 1-D product between them.
We thus explored the precision of this approximation in different scenarios (Fig. 2). To do so, we randomly generated groups of nodes of different order with different values of {N i } (out of N total hyperlinks) and different values for the number of collective hyperlinks N 1,2,...,n . Specifically, for each size n we extracted n values of N i from the uniform distribution U (1, f N ), where f is the fraction of maximum number of hyperlinks out of the total N that we allow a node to appear in. We then extracted N 1,2,...,n from the uniform distribution U (1, min(N i )). We then used the values (N 1,2,...,n , N i , N ) to extract a p-value with the exact approach of Eq. 3 and the approximated one of Eq. 4. We varied the size of the simplices to validate in the interval [1,4]. Each column of Fig. 2 presents results with a different size. We also varied the total number of interactions  Fig. 2 presents results with a different value of f . We find that across the different settings, the approximated p-value tends to overestimate the effective probability that the number of interactions within a group is compatible with our null hypothesis. This can lead to False Negative, i.e. to miss groups that are significant but for which we overestimate the actual p-value. However, this overestimation is pronounced only for values of f close to 1, which means cases in which a system is made of a small amount of nodes each of which is responsible, on average, of a significant portion of the complete number of hyperlinks. In cases in which a single node is at maximum present in a smaller fraction of the overall interactions, the approximation works well and the p-values computed using the two approaches fall on the diagonal. In all the applications that we are going to showcase in the following, we have used the approximated p-value after verifying that the maximum value of f for each node in all investigated systems was below 0.1.

A. Definition and investigation of benchmarks
To test the effectiveness of our approach, we create benchmark hypergraphs made of N nodes in which we implant T simplices for each size n. In this work, we explored N = {100, 200, 500} and T = {100, 200, 400}, with n ranging in the interval [2,4]. For each simplex, we create l hyperlinks, where l is sampled from a Binomial distribution with p bin = 0.5 and n bin = 20. Each hyperlink contains the n nodes of the simplex plus a randomly selected set of the other N − n nodes whose exact size is sampled from the uniform distribution U (0, n max ). In our simulations, we set n max = 6. In this way the n nodes of a simplex interact l times in different hyperlinks, which always contain the n core nodes plus a set of randomly selected ones. Moreover, we introduce a pa- rameter c that controls the level of simplicial closure at each size. Specifically, c is the fraction of groups of size n that constitute a clique of size n + 1. For n = 2, c corresponds to the concept of triadic closure and represents the fraction of dyads that are closing a triangle. For each realization of our benchmark, we use (i) the pairwise method of Statistically Validated Networks (SVN) [39], (ii) the method of Statistically Validated Hypergraphs (SVH), that only validates observed hyperlinks of fixed size [48], and (iii) our new method of Statistically Validated Simplices (SVS) to extract the simplices rejecting the null hypothesis. To validate our results, we compute the True Positive Rate (TPR), that quantifies the ratio of true groups that we are able to identify, and the False Discovery Rate (FDR) that quantifies the ratio of erroneously detected groups that the different methods select. Specifically, TPR= T P T P +F N and F DR = F P F P +T P , where TP are the True Positives (in our case, simplices correctly identified), FN are the False Negatives (significant simplices present in the system but not identified by our method) and FP are the False Positives (simplices wrongly identified as significant). Fig. 3a shows TPR as a function of c for N=200 and T=100, but the observed behavior is independent of the specific values of N and T. Each point is the median of 100 benchmark realizations. The color bands represent the interval between the 10 t h and the 90 t h percentile. The performance of SVS and SVH is not affected by c, as both are methods designed to take into account the higher-order nature of interactions. In fact, if 3 nodes i, j, k constitute a clique of (significant) pairwise interactions but never interact altogether, SVS and SVH will not statistically validate the corresponding simplex because N ijk = 0 even if N ij , N ik and N jk are not. As a consequence, SVS and SVH will not miss {ij}, {jk} and {ik} as significant simplices. Conversely, the performance of SVN significantly decreases with c specifically because the SVN puts together nodes that never interacted altogether and misses the smaller true simplices. In our example, SVN will erroneously identify {ijk} as a significant simplex and will miss that the significant ones are the couples {ij}, {jk} and {ik}. We also observe that SVS significantly outperforms SVH. This happens because in our benchmark hypergraphs we allow the n nodes of a simplex to interact in larger groups. While SVS counts all the hyperlinks in which the n nodes appear, the SVH approach misses all those in which there are additional, randomly selected nodes. For this reason, SVH is more likely to underestimate the significance of a simplex and thus have a lower value of TPR. Fig. 3b plots the FDR as a function of c. SVH and SVS almost never detect wrong groups as they distinguish between interactions of different order. SVN instead has a significantly higher value of FDR, which worsens for high c because the number of simplices erroneously aggregated at larger sizes increases, as pointed out in our example discussed previously. To better characterize the difference in performance between SVS and SVH, we add to the benchmark a parameter f that regulates the fraction of times in which the n nodes that constitute a simplex interact in hyperlinks of size n > n. The larger f , the more we expect SVH to underestimate the number N 12...n of collective interactions and thus will likely fail to validate the true simplex. In Fig. 3c we show TPR as a function of f . Indeed, we see that SVH significantly worsens its performance when f increases. Fig. 3d shows the behavior of FDR as a function of f . SVH does not detect a significantly larger ratio of false simplices until f is so large that it detects very few positives (compare with TPR for f ≥ 0.8, all of which are false.
B. Application to real data: Baskets of co-purchased goods at Walmart We analyse a dataset that tracks 69,906 co-purchases of 88,860 products at Walmart [64]. This system can be mapped into a hypergraph where hyperlinks represent copurchases. In this context, we use the SVS approach to extract simplices of products that are inherently bought together. When extracting the Statistically Validated Simplices we exclude from the hypergraph all baskets of size larger of 10. Indeed, each large basket of size k which is not statistically significant produces a large number of candidate simplices of smaller size k -specifically it produces k k candidates. Given that larger baskets have small number of co-occurrence, we cut the system at size 10 in order to avoid the explosion in the number of hyperlinks with low chance of getting validated that we check. Fig. 4a shows the number of validated (blue) and not validated (orange) simplices for each size. The number of validated simplices exponentially decays with size, and, for large n, it is several orders of magnitude smaller than the simplices that are observed at least once but are not rejecting the null hypothesis. The dataset also reports the category of each product. By using this information, we define a proxy of "purity" of co-purchased goods. Purity is a measure that goes from 0 (if all products in a group belong to different categories) to 1 (if they all belong to the same category), and is defined as P = n n−1 (1 − C n ), where n is the size of the simplex and C is the number of unique categories to which the n products belong. Fig. 4b shows the "purity" with respect to category classifications of the validated (blue) and not validated (orange) simplices as a function of their size. In the figure, the solid lines plot the mean and the transparent bands the error of the mean. We find that the SVS of any size are more homogeneous with respect to category than the ones that are not validated.

C. Co-traded stocks belonging to the S&P100 index
We investigate a dataset tracking tick-by-tick transactions performed by 96 stocks that belong to the S&P100 index and are traded on any US venue in the month 12/2020 [65]. The data comes with a temporal resolution of 1 ms. By applying our method to this dataset, we detect SVSs of stocks that are co-traded in the investigated period in excess to what expected by a null model taking into account the heterogeneity observed in stock trading. Also here, before extracting the SVSs we exclude from the hypergraph all simplices larger than 10 stocks. We perform our analysis at a daily scale and then aggregate the results obtained for different days. Fig. 5a shows the number of validated (blue) and not validated (orange) simplices for each size. Both for validated and not validated simplices, Fig. 5b shows the mean purity with respect to the economic sector of stocks (solid lines) and its error (transparent bands) as a function of size. We find that validated simplices (blue line) have higher purity at smaller size but the difference disappear in larger ones. Fig. 5c shows the mean correlation (solid line) and its error (transparent bands) between couples of stocks in validated (blue) and not validated (orange) simplices. Correlation is computed on time series of financial returns extracted at a resolution of 1 minute. We find that stocks belonging to validated simplices have higher correlation patterns than those that are not validated.

III. DISCUSSION
The recent surge in interest towards higher-order networks, together with the increased availability of large amount of data on several real world systems, have lead to the introduction of the first filtering tools for hypergraphs [45,48]. However, the task of detecting statistically validated units of collectively interacting agents is multi-faceted, and the multidimensional nature of higherorder interactions calls for the design of different validation approaches. The existing methods have focused on the extension of validating approaches from pairwise in- teractions to hyperlinks [45] and the detection of statistically validated groupwise interactions that occur only at a specific size [48]. Here, we introduced a different and more general framework for the detection of Statistically Validated Simplices, which represent sets of nodes that appear together in interactions of any size. The computational efforts needed for the statistical validation are demanding but accessible with ordinary PCs for groups up to a size of about 10 elements. To speed up the computation we have proposed an approximated procedure that work well in a large number of cases as shown in our investigation of benchmarks. Using properly designed hypergraph benchmarks, we show that our approach is highly effective on higher-order systems (i) in which informative interactions among n core nodes are likely to be diluted into larger hyperlinks that feature occasional nodes and (ii) with a high degree of simplicial closure. We also show how, on real data, our approach allows to detect simplices whose nodes are characterized by high levels of similarity, which might help in shredding light on the generative process of the underlying hypergraphs.

Code Availability
The code to extract Statistically Validated Simplices is available at https://github.com/musci8/SVH.