Network resampling for estimating uncertainty

Next, we describe in detail the three subsampling procedures we study. While there are other subsampling strategies and some may have advantages in specific applications, we focus on these three because they are popular and general.

Before we proceed to the data-driven method of choosing the settings for resampling, we summarize the algorithm we employ for constructing confidence intervals for network statistics based on a given resampling method with a given sampling fraction $q$. This is a general bootstrap-style algorithm that can be applied to any network statistic $T$ computable on a subgraph.

For any method of network subsampling, we have to choose the sampling fraction; for consistency across different methods, we focus on the choice of the node pair fraction $q$, which can be converted to $p$ for node or row sampling. In the classical bootstrap analogue, the $m$-out-of-$n$ bootstrap, consistency results have been obtained for $m\rightarrow\infty$ and $m/n\rightarrow 0$ as $n\rightarrow\infty$. These do not offer us much practical guidance for how to choose $m$ for a given finite sample size $n$. The situation is similar for network bootstrap, with no practical rules for choosing $q$ implied by consistency analysis.

The iterated bootstrap has been commonly used to calibrate bootstrap results since it was proposed and developed in the 1980s [16, 1, 17]. An especially relevant algorithm for our purposes is the double bootstrap approach, originally proposed to improve coverage of bootstrap intervals [28]. A similar double bootstrap approach was proposed by [22] for $m$-out-of-$n$ bootstrap, to correct a discrepancy between the nominal and the actual coverage for bootstrap confidence intervals due to a different sample size $m$.

A generic double bootstrap algorithm proceeds as follows. Given the original i.i.d. sample $\mathcal{X}=\{X_{1},\dots,X_{n}\}$ and an estimator $\hat{\theta}$ based on $\mathcal{X}$, consider a size $n$ bootstrap sample $\mathcal{X}^{*}$, drawn from $\mathcal{X}$ with replacement, and the corresponding estimator $\hat{\theta}^{*}$. If the distribution of $\hat{\theta}$ can be approximated by the distribution of $\hat{\theta}^{*}$ conditional on $\mathcal{X}$, a confidence interval of nominal level $1-\alpha$ can be constructed as $I(\alpha)=[\hat{t}_{\alpha/2},\hat{t}_{1-(\alpha/2)}]$, with $\hat{t}_{\alpha}$ defined as $\sup\{t:\mathbf{P}(\hat{\theta}^{*}\leq t\mid\mathcal{X})\leq\alpha\}$. This confidence interval typically has correct coverage up to $O(n^{-1})$, and double bootstrap proposes to improve this bound by using $I(\beta_{\alpha})$ instead, where $\beta_{\alpha}$ is defined as the solution to $\mathbf{P}(\theta\in I(\beta_{\alpha}\mid\mathcal{X}))=1-\alpha$ and can be approximated by $\hat{\beta}_{\alpha}$, the solution to $\mathbf{P}(\hat{\theta}\in I(\hat{\beta}_{\alpha}\mid\mathcal{X}^{*})\mid\mathcal{X})=1-\alpha$. In practice, this is implemented by sampling multiple size $n$ $\mathcal{X}^{**}$ draws from each $\mathcal{X}^{*}$ with replacement, and estimating $\mathbf{P}(\hat{\theta}\in I(\beta_{j}\mid\mathcal{X}^{*})\mid\mathcal{X})$ as the proportion of confidence intervals based on $(1-\beta_{j})$ percentiles of $\mathcal{X}^{**}$ that cover $\hat{\theta}$ for a set of candidate $\beta_{j}$’s. We then choose $\hat{\beta}_{\alpha}=\beta_{j}$ such that $\mathbf{P}(\hat{\theta}\in I(\beta_{j}\mid\mathcal{X}^{*})\mid\mathcal{X})$ is the closest to $1-\alpha$ among all candidate $\beta_{j}$’s.

Since different choice of subsampling fraction $q$ will lead to different level to coverage, here we propose an analogous double bootstrap method for network resampling. Compared to the i.i.d. case, instead of using size $n$ resamples sampled with replacement, we will generate resamples using network subsampling methods described above and instead of choosing an appropriate $\alpha$, we aim to choose an appropriate $q$. The algorithm is detailed as follows:

Input: a graph $G$, a network statistic $T(G)$, a subsampling method, a set of candidate fractions $q_{1},\dots,q_{J}$, and the number of bootstrap samples $B$.

For $j=1,\dots,J$,

1. 1.

   For $b=1,\dots,B$,

   1. (a)

      Split all nodes randomly into two sets $\mathcal{V}_{1}^{b}$, $\mathcal{V}_{2}^{b}$.

   2. (b)

      On the subgraph induced by $\mathcal{V}_{1}^{b}$, calculate the estimate $\hat{T}_{b,j}^{(1)}$.

   3. (c)

      On the subgraph induced by $\mathcal{V}_{2}^{b}$, run Algorithm 1 with $q=q_{j}$ to construct a confidence interval $R_{b,j}$ for $T_{b,j}^{(2)}$

2. 2.

   Calculate empirical coverage rate $\pi_{j}=\frac{1}{B}\sum_{b=1}^{B}\mathbf{1}\{\hat{T}_{b,j}\in R_{b,j}\}$.

Output: $\hat{j}=\arg\min(\pi_{j}-(1-\alpha))^{2}$.

One key difference from the i.i.d. case is that we estimate empirical coverage by checking if the confidence interval constructed on each subsample covers $\hat{T}_{b,j}^{(1)}$, not $\hat{T}(G)$, the estimate on the original observed graph $G$. This is because with network subsampling, we cannot preserve the original sample size $n$, and if we are going to compare network statistics on two subsamples, we need to match the graph sizes (in this case $n/2$).

Count statistics are functions of counts of subgraphs in network, and are commonly viewed as key statistical summaries of a network [3, 4], analogous to moments for i.i.d.  data. Here we will examine a simple and commonly studied statistic, the density of triangles. A triangle subgraph is three vertices all connected to each other. Triangle density is frequently studied because it is related to the concept of transitivity, a key property of social networks where two nodes that have a “friend” in common are more likely to be connected themselves. We will apply resampling to construct confidence intervals for the normalized triangle density for a given network.

We generate the networks from the SBM models with number of communities $K$ set to either 1 (no communities, the Erdös-Rényi graph) or 3, with equal community sizes. We consider two values of $n$, 300 and 600. The matrix of probabilities $B$ for the SBM is defined by, for all $k,l=1,\dots,K$,

The difficulty of community detection in this setting is controlled by the edge density $\rho$ and the ratio $t=\gamma_{1}/\gamma_{2}$. We set $t=5$, and vary the expected edge density $\rho$ from 0.01 to 0.1, with higher $\rho$ corresponding to more information available for community detection.

For a fair comparison, we match the proportion of adjacency matrix observed, $q$, in all subsampling schemes, and vary $q$ from 0.2 to 0.8. Figures 2 and 3 show the widths and coverage rates for the confidence intervals constructed by the three subsampling procedures over the full range of $q$, as well as the value of $q$ chosen by our proposed double bootstrap procedure.

There are several general conclusions we can draw from the results in Figures 2 and 3.

(1) The confidence intervals generally get shorter as $q$ increases, since the subsampled graphs overlap more, and this dependence reduces variance. The coverage rate correspondingly goes down as $q$ goes up. Still, all resampling procedures achieve nominal or higher coverage for $q\leq 0.5$. The one exception to this is node pair sampling at the lowest value of $q$ with the lowest edge density, where most subsampled graphs are too sparse to have any triangles. In this case the procedure fails, producing estimated zero triangle density and confidence intervals with length and coverage both equal to zero.

(2) For a given value of $q$, a larger effective sample size, either from higher density $\rho$ or from a larger network size $n$, leads to shorter intervals, as one would expect.

(3) Node sampling produces shorter intervals than the other two schemes, while maintaining similar coverage rates for smaller values of $q$. The coverage rate falls faster for node sampling for larger values of $q$ outside the optimal range.

(4) The double bootstrap procedure selects a value of $q$ that achieves a good balance of coverage rate and confidence interval width, maintaining close to nominal coverage without being overly conservative. It does get misled by the pathological sparse case with no triangles, highlighting the need of not deploying resampling algorithms without checking that there is enough data to preserve the feature of interest (triangles in this case) in the subsampled graphs.

Figure 4 shows the value of $q$ chosen by the double bootstrap procedure as a function of graph size and edge density; both appear to not have much effect on the value of $q$ within the range considered. This suggests that as long as the structure of a network is preserved, a reasonable value of $q$ could be chosen on a smaller subsampled network, resulting in computational savings with little change in accuracy. And we can see that the $q$ chosen for node sampling is always smaller than that for the other two subsampling methods, so node sampling can be more computationally efficient, giving smaller subsampled graphs.

Finding communities in a network often helps understand and interpret the network patterns, and the number of communities $K$ is required as input for most community detection algorithms [2, 12]. A number of methods have been proposed to estimate $K$: as a function of the observed network’s spectrum [21], through sequential hypothesis testing [8], or subsampling [25, 9]. Subsampling has been shown to be one of the most reliable methods to estimate $K$, but little attention has been paid to choosing the subsampling fraction and understanding its effects. Here we investigate the performance of the three subsampling schemes on this task and the choice of the subsampling fraction on the task of understanding uncertainty in $K$, beyond the earlier focus on point estimation. Again, we match the proportion of adjacency matrix observed, $q$, across all subsampling schemes.

For this task, we generate networks from an SBM with $K=3$ equal-sized communities. We vary the number of nodes $n$, the expected edge density $\rho$, and the ratio of within and between communities edge probabilities $t$, which together control the problem difficulty.

For subgraphs generated by node sampling and row sampling, we estimate the number of communities using the spectral method based on the Bethe-Hessian matrix [21, 33]. For a discrete variable like $K$, confidence intervals make little sense, and thus we look at their bootstrap distributions instead. These distributions are not symmetric, and rarely overestimate the number of communities. For subgraphs generated from node pair sampling, there is no complete adjacency matrix to compute the Bethe-Hessian for, and computing it with missing values replaced by zeros is clearly introducing additional noise. We instead estimate the number of communities from node-pair samples using the network cross-validation approach [25], first applying low-rank matrix completion with values of rank $K$ ranging from 1 to 6, and choosing the best $K$ by maximizing the AUC score for reconstructing the missing edges.

Figure 5 and 6 shows the results under different subsampling schemes. Since $K$ is discrete, instead of considering coverage rate, we compare methods by the proportion of correctly estimated $K$. It is clear that whenever the community estimation problem itself is hard (the lowest values of $\rho$, $t$, and $n$), the number of communities is also hard to estimate. In these hardest cases, $K$ tends to be severely underestimated by node and row sampling, and varies highly for node pair sampling, and there is no good choice of $q$. In all other scenarios, our proposed algorithm does choose a good value of $q$ balancing the bias and variability of $\hat{K}$ by maintaining the empirical coverage rate above 90% and provide a set of candidate rank estimates.

In many network applications, we often observe covariates associated with nodes. For example, in a school friendship network survey, one may also have the students’ demographic information, grades, family background, and so on [18]. In such settings, the network itself is often used as an aid in answering the main question of interest, rather than the primary object of analysis. For instance, a network cohesion penalty was proposed as a tool to improve linear regression when observations are connected by a network [24]. They consider the predictive model

where $X\in\mathbb{R}^{n\times p}$, $Y\in\mathbb{R}^{n}$, and $\beta\in\mathbb{R}^{p}$ are the usual matrix of predictors, response vector, and regression coefficients, respectively, and $\mathbf{\alpha}\in\mathbb{R}^{n}$ is the individual effect for each node associated with the network. They proposed fitting this model by minimizing the objective function

where $\lambda_{1}$ is a tuning parameter, and $L=D-A$ is the Laplacian of the network, with $D=\text{diag}(d_{1},\dots,d_{n})$ the degree matrix. The Laplacian-based network cohesion penalty has the effect of shrinking individual effects of connected nodes towards each other, ensuring more similar behavior for neighbors, and improving prediction [24].

This method provides a point estimate of the coefficients $\beta$, but no measure of uncertainty, which makes interpretation difficult. Our subsampling methods can add a measure of uncertainty to point estimates $\hat{\beta}$. Since different subsampling methods will lead to different numbers of remaining nodes and thus different sample sizes for regression, we only replace the penalty term in (1) with $\lambda_{1}\alpha^{T}\tilde{L}\alpha$, where $\tilde{L}$ is the corresponding Laplacian of the subsampled graph, and retain the full size $n$ sample of $Y$ and $X$ for regression, to isolate the uncertainly associated with the network penalty. Alternatively, one could resample $(X,Y)$ pairs together with the network to assess uncertainty overall. If we view the observed graph $A$ as a realization from a distribution $F$, then the subsampled graphs under different subsampling methods can also be considered as a smaller sample from the same distribution $F$ and the estimation based on these smaller samples should still be consistent.

In addition to the three subsampling methods, we also include a naive bootstrap method. Let $A$ be the original adjacency matrix with $n$ nodes. Sample $n$ nodes with replacement, obtaining $n_{1},n_{2},\dots,n_{n}$, and create a new adjacency matrix $\tilde{A}$ by setting $\tilde{A}_{ij}=1$ if $A_{n_{i}n_{j}}=1$ or $n_{i}=n_{j}$, otherwise, $\tilde{A}_{ij}=0$. Note that with this method we cannot isolate the network as a source of uncertainty, and thus the uncertainty assessment will include the usual bootstrap sample variability of regression coefficients.

We generate the adjacency matrix $A$ as an SBM with with $K=3$ communities, with 200 nodes each. The individual effect $\alpha$ is generated from the normal distribution $N(c_{k},\sigma_{\alpha}^{2})$, where $c_{k}$ are $-1$, $0$, and $1$ for nodes from communities $1$, $2$, and $3$, respectively. We vary the ratio of within and between edge probabilities $t$ and the standard deviation $\sigma_{\alpha}$, to obtain different signal-to-noise ratios. Regression coefficients $\beta_{j}$’s are drawn from $N(1,1)$ independently, and each response $y_{i}$ is drawn from $N(\alpha_{i}+\beta^{T}x_{i},1)$. Figures 7, 8 and 9 show the results from this simulation.

Once again, as $q$ increases, the confidence intervals get narrower as subgraphs become more similar to the full network. However, for node and row sampling the coverage is also lower for the smallest value of $q=0.1$, since the variance of bootstrap estimates of $\beta$ is larger for smaller $q$, and the center of the constructed confidence interval gets further from the true $\beta$. This may be caused by the fact that only a small fraction of nodes contribute to the penalty term. Node pair sampling is different: it gives better performance in terms of mean squared error for nearly all choices of $q$, but since its confidence intervals are much narrower than those of node and row sampling, their coverage rate is very poor. All subsampling methods are anti-conservative and can only exceed the nominal coverage rate with small $q$. One possible explanation is that in a regression setting, the impact of node covariates dominates any potential impact of the network. In all cases, Algorithm 2 chooses $q=0.1$, which is either the only choice of $q$ that exceeds the nominal coverage rate or the choice of $q$ that gives a coverage rate close to the nominal coverage rate.

We now consider the case when the dimension of node covariates is larger. The LASSO in [36] is a popular technique used for high dimensional i.i.d. data. By minimizing

the resulted estimator $\hat{\beta}$ will have some entries shrunk to zero, and the set $S=\{i:\hat{\beta}_{i}\neq 0\}$ is known as the activated set. However, we need to balance the goodness-of-fit and complexity of the model properly by tuning $\lambda_{2}$, and model selection using cross validation is not always stable. In [29], the authors proposed stability selection, which uses subsampling to provide variable selection that is less sensitive to the choice of $\lambda$.

Here we propose a similar procedure for network data. With a graph $G$ observed, we first generate $B$ copies of $g_{b}$ using one of the resampling schemes. Then we minimize the objective function

over $g_{b}$ to construct activated set $S_{b}$, and for each $1\leq j\leq p$, calculate $\frac{1}{B}\sum\mathbf{1}\{j\in S_{b}\}$, the selection probability of each predictor.

The simulation setup is going as follows: A graph with 3 blocks with 200 nodes in each community is generated through SBM. Again we vary the ratio of within and between communities edge probability and the edge density of the graph. The individual effect $\alpha$ is generated from $N(c_{k},\sigma_{\alpha})$, where $c_{k}$ are -1, 0, 1 for nodes from different communities respectively. 25 coefficients $\beta_{j}$ are drawn from $N(1,1)$ independently, while the rest 75 $\beta_{j}$ are set to 0. Then $y_{i}$ is drawn from $N(\alpha_{i}+\beta^{T}x_{i},1)$. The simulation result is shown in Figure 10, 11, 12.

The three resampling schemes give different performances in AUC when we vary the subsampling proportion $q$: the performance of $p$ sampling is not sensitive to the choice of $q$, while the performances of row sampling and node pair sampling improve or worsen respectively as $q$ increases.

Similar to what we observed in the last section, although the subsampling step does bring in uncertainty into the inference, in a regression setting, node covariates often play a more important role and the performance of resampling only the network structure is not as good as expected.