1 Introduction↩︎

Link prediction [1]–[4] is a fundamental problem in network data analysis, with applications in sociology, econometrics, and biology, where real-world networks are often only partially observed. While early work focused on accurately predicting individual links, the growing demand for simultaneous prediction of multiple links has brought concerns about reproducibility and statistical validity to the forefront. Framing this task as a multiple hypothesis testing problem has led to the adoption of false discovery rate (FDR) control as a key tool. However, FDR control for link prediction in networks presents unique challenges owing to complex edge dependencies and unknown, heterogeneous patterns of missingness. To address these challenges, which remain inadequately resolved in the literature, we introduce conformal link prediction (clp), a new method for simultaneous link prediction with provably FDR control.

Let \({\mathbf{A}}\in \mathbb{R}^{n\times n}\) denote the adjacency matrix of a directed network with \(n\) nodes and no self-loops, where entries \(A_{i,j}\) and \(A_{j,i}\) are not necessarily equal. Missingness in \({\mathbf{A}}\) is described by a binary matrix \({\mathbf{M}}\in \{0,1\}^{n\times n}\), where its entry \(M_{i,j}=1\) if \(A_{i,j}\) is missing, and \({\mathbf{M}}\) can be any deterministic or random binary matrix. We define \({\mathbf{M}}^c=\boldsymbol{1}_{n\times 1}\boldsymbol{1}_{n\times 1}^{{ \mathrm{\scriptscriptstyle T} }}-{\mathbf{M}}\). Throughout, we only adopt the following standard assumption.

Assumption 1 is closely related to the classical missing-completely-at-random (MCAR) condition in the missing data literature [6]–[8]. In our graphon setting, it requires the missingness pattern \({\mathbf{M}}\) to be independent of the realized network \({\mathbf{A}}\), while still allowing heterogeneous or structured missing patterns. This differs from missing-at-random (MAR) mechanisms, where missingness may depend on the latent node variables \(\{\xi_i\}_{i=1}^n\) but not directly on the realized link values conditional on these latent variables. It also excludes missing-not-at-random (MNAR) mechanisms, where missingness may depend directly on the unobserved value \(A_{i,j}\).

We define the test set as the collection of unobserved links: \({\cal D}_{\rm test}:=\{(i,j): i,j \in [n]=\{1,\ldots, n\}, M_{i,j}=1\}\), where \({\cal D}_{\rm test}\) is random if \({\mathbf{M}}\) is random. For simplicity of presentation, we focus on predicting all links in \({\cal D}_{\rm test}\), while our method can be applied to any subset. Following [9], we formulate link prediction as a multiple testing problem with random hypotheses: \[H_{i,j}: A_{i,j} \leq c_{i,j}~ \rm for(i,j) \in {\cal D}_{\rm test}, \label{def::hypothesis-testing}\tag{1}\] where \(c_{i,j}\)’s are user-specified. For example, in the international trade flow network [10], if the goal is to determine whether any trading activity exists between two countries, \(c_{i,j}\) can be set to zero. Alternatively, if the focus is on detecting higher volumes of bilateral trade, \(c_{i,j}\) can be set to an appropriate positive value.

For the rejection set \(\mathcal{R}:=\{{(i,j) \in {\cal D}_{\rm test}}: \textrm{Reject }H_{i,j}\}\), the false discovery proportion (FDP) is defined as \(\textrm{FDP}= (|\mathcal{R}| \vee 1)^{-1} \sum_{(i,j)\in \mathcal{R}} \mathbb{1}_{[A_{i,j}\leq c_{i,j}]}\), based on given \(\{c_{i,j}\}_{(i,j) \in {\cal D}_{\rm test}}\) and observed data \(\big\{{\mathbf{A}}\circ {\mathbf{M}}^c, {\mathbf{M}}\big\}\), where “\(\circ\)” denotes the element-wise product. Our objective is to control the conditional FDR given \({\mathbf{M}}\) at a pre-specified significance level, and accordingly treating the test set as fixed. Clearly, the marginal FDR control is thus also valid since one can marginalize the conditional FDR control guarantee over \({\mathbf{M}}\). Specifically, we define \(\textrm{FDR}=\mathbb{E}(\textrm{FDP}|{\mathbf{M}})\). A major challenge for FDR control in classical link prediction methods is that the unknown and potentially non-uniform missing patters lead to non-exchangeable observations \({\mathbf{A}}\circ {\mathbf{M}}^c\), which hinders uncertainty quantification. Moreover, complex dependencies among predicted links further render most traditional FDR control procedures inapplicable.

1.1 Related literature↩︎

Link prediction has been extensively studied over recent decades [11]–[13]. An alternative perspective treats the network as a random matrix or a realization of an unobservable graphon, leading to approaches based on matrix estimation or completion [14]–[16] and graphon estimation [17], [18]. However, most existing studies either lack uncertainty quantification or rely on strong modeling assumptions. To achieve our goal of distribution-free uncertainty quantification, we bridge two key ingredients. The first is conformal prediction [19], a powerful model-agnostic framework that enables flexible applications, including recent adaptations to multiple testing via the construction of conformal \(p\)-values [9], [20], [21]. The second ingredient is the e-BH procedure [22], which ensures valid FDR control under complex dependence among conformal \(p\)-values. Based on e-values [23], this procedure has proven effective in handling complex dependence among statistics, with recent successful applications in derandomized knockoffs [24] and e-value-based conformal prediction [25], [26].

Conformal prediction for matrix data was pioneered by [27] and [5]. [27] considered random missingness and is therefore more closely related to our setting, whereas [5] examined arbitrary and potentially adversarial missing patterns. However, neither study addresses the task of predicting multiple links. Recently, [28] introduced conformal link prediction with FDR control, but the method is limited to unweighted networks and assumes a “double-standard" sampling scheme, requiring uniform missing rates across all links. Under this assumption, subsequent work by [29] sought to provide theoretical justification for [28]. This uniformity condition effectively imposes the exchangeability on all links with values exceeding \(c_{i,j}\), simplifying the analysis but limiting practical applicability. [30] proposed a structured splitting scheme for simultaneous prediction of a small number of links, without considering the multiple testing. [31] construct prediction intervals for individual links with marginal coverage guarantees, rather than controlling the FDR for multiple dependent links simultaneously. To date, valid FDR control for multiple link prediction under general network settings remains an open challenge.

1.2 Our contributions↩︎

We develop a new conformal inference framework for multiple network link prediction with provable FDR control. Unlike prior work, our method is tailored to the complex structure of partially observed networks, characterized by heterogeneous missingness and intricate dependencies among predictions. The approach combines novel conformal calibration schemes with a principled global aggregation strategy. Our main contributions are threefold.

First, we introduce a conformal inference method tailored to row-wise link prediction. A key innovation is a carefully constructed multi-splitting scheme that partitions each row’s test set and adaptively allocates calibration sets. This enables the construction of exchangeable nonconformity scores and valid conformal \(p\)-values, even under unknown and nonuniform missingness. The method achieves FDR control for any node sizes \(n\) through a nontrivial adaptation of the Benjamini–Hochberg (BH) procedure [32], applied in a row-wise fashion.

Second, we propose a novel aggregation technique that transforms the row-wise test decisions, obtained via the BH procedure, into valid e-values, and applies the e-BH procedure for FDR control. This step requires careful design to accommodate unknown yet arbitrary dependence across predicted network links, a challenge not addressed by standard conformal or e-value methods. We further introduce a derandomization step to enhance stability, along with an e-value inflation strategy to improve power, both tailored to the unique demands for predicting multiple network links.

Third, our method is broadly applicable: it extends naturally to unweighted, undirected, and bipartite networks, requires no knowledge of the missingness mechanism, and provides finite-sample guarantees. Our framework advances the scope of conformal inference in network analysis, extending beyond existing approaches such as [33] and [27]. Numerical studies confirm the practical effectiveness of the proposed method.

2 Conformal link prediction with FDR control↩︎

2.1 Link prediction under unknown missing probabilities↩︎

We adopt the general weighted graphon model [34] as the base for our method for its flexibility and generality. We assume entry \(A_{i,j}\) of the adjacency matrix \({\mathbf{A}}\) is drawn from a distribution parameterized by \(f(\xi_i,\xi_j)\), such as Bernoulli and Gaussian. That is, \[\begin{align} A_{i,j} \sim f(\xi_i,\xi_j), \rm ~for~i\neq j \in [n], ~~ \rm where \{\xi_i\}_{1\leq i\leq n} \stackrel{\rm i.i.d.} \sim \mathrm{Uniform}[0,1]. \label{def::exch-model} \end{align}\tag{2}\] Intuitively, model 2 first samples nodes independently from a hyper-population of “all kinds of nodes", then generates links according to the function \(f\). This framework encompasses a broad class of network models, such as in the social relations model [35], the conditionally independent dyad model [36], and the random-effects model [37].

Under model 2 , for any permutation \(\pi(\cdot)\) of node labels, we have \[\{A_{i,j}\}_{i\neq j \in [n]} \stackrel{d}{=} \{A_{\pi(i),\pi(j)}\}_{i\neq j \in [n]}, \label{def::exch-property}\tag{3}\] which is known as the exchangeability property of \({\mathbf{A}}\). Importantly, the exchangeability in 3 pertains to the not-fully-observable complete \({\mathbf{A}}\), not the observed \({\mathbf{A}}\circ {\mathbf{M}}^c\). Note that our notion of exchangeability in 3 differs from that in prior work. Previous studies [27], [33] rely on exchangeability induced by the random \({\mathbf{M}}\) where they assume \(M_{i,j}\sim\operatorname{Ber}(q_{i,j})\), and therefore require either uniform missing rates \({q_{i,j}}\) or their accurate estimations under the “weighted exchangeability" framework. In contrast, we model \({\mathbf{A}}\) as a random graph with inherent exchangeability as defined in 3 . With the mild additional assumption of exchangeability of the complete \({\mathbf{A}}\), we waive the need for prior knowledge of the potentially non-uniform missing rates \({q_{i,j}}\)’s, a key restrictive assumption in earlier works.

To address the challenge of unknown missing patterns, we propose a divide-and-conquer approach to break the problem into sub-tasks. First, for each row \(i_0\), we construct “local tests” where exchangeability holds within a row-wise subset of observations, which leads to valid “local” FDR controls. Then, the local test decisions are aggregated into e-values to construct an aggregated test of the hypotheses in 1 for all unobserved links, with valid FDR control.

2.2 Step 1: breaking down original test into smaller local tests↩︎

Our first step is to design a local test, which focuses only on row \(i_0\), by constructing exchangeable conformal \(p\)-values [9]. Let \({\cal I}_{\rm test}^{i_0}=\{(i_0,j): j\in [n], M_{i_0,j}=1\}\) denote the column indices of missing entries in row \(i_0\), corresponding to missing links for node \(i_0\). Then, the corresponding set of local hypotheses is \[H_{i_0,j}: A_{i_0,j} \leq c_{i_0,j}, ~~ \rm for j \in {\cal I}_{\rm test}^{i_0}. \label{def::hypothesis-testing-local}\tag{4}\] We randomly split the observed entries in row \(i_0\), that is, \(\{A_{i_0,j} : j \notin {\cal I}_{\rm test}^{i_0}\}\), into a training set \(\{A_{i_0,j} : j \in {\cal I}_{\rm train}^{i_0}\}\) and a calibration set \(\{A_{i_0,j} : j \in {\cal I}_{\rm calib}^{i_0}\}\), as illustrated in Figure 1.

Under the random missing mechanism in Assumption 1, the observed entries in the same row can be viewed as a randomly selected subset of exchangeable random variables. It turns out that, as a consequence of the exchangeability of latent factors \(\{\xi_i\}_{1\leq i\leq n}\), the entries in \(\{A_{i_0,j}\}_{j \in {\cal I}_{\rm calib}^{i_0}\cup {\cal I}_{\rm test}^{i_0}}\) are exchangeable, as established in the following proposition.

For each \(j \in {\cal I}_{\rm calib}^{i_0}\cup {\cal I}_{\rm test}^{i_0}\), we select data from observed links using a carefully designed strategy, detailed later in Section 2.3, to construct a link prediction \(\widehat A_{i_0,j}\) such that \[\begin{align} \{\widehat A_{i_0,j}\} \rm are exchangeable forj\in {\cal I}_{\rm calib}^{i_0}\cup{\cal I}_{\rm test}^{i_0}. \label{eq:est-cond} \end{align}\tag{5}\] Based on 5 , we construct the nonconformity scores [9] as follows \(S(\widehat A_{i_0,j}; z) := z - \widehat A_{i_0,j}\) for \(j \in {\cal I}_{\rm calib}^{i_0}\cup {\cal I}_{\rm test}^{i_0}\), where where \(S(\widehat A_{i_0,j}; z)\) is monotone with respect to the second argument by definition. Intuitively, \(S(\widehat A_{i_0,j}; z)\) quantifies how well the value \(z\) conforms to the predictions. Using the nonconformity scores, for each \(j_1 \in {\cal I}_{\rm test}^{i_0}\), we define \(p\)-values for entries \(\{(i_0,j_1):j_1 \in {\cal I}_{\rm test}^{i_0}\}\) as \(\widehat{\mathfrak p}_{i_0,j_1}= (1 + |{\cal I}_{\rm calib}^{i_0}|)^{-1}[1 +\sum_{j' \in {\cal I}_{\rm calib}^{i_0}}\mathbb{1}\{ S(\widehat A_{i_0,j'}; A_{i_0,j'}) < S(\widehat A_{i_0,j_1}; c_{i_0,j_1}) \}]\). We then apply the BH procedure to \(\big\{\widehat{\mathfrak p}_{i_0,j}: j \in {\cal I}_{\rm test}^{i_0}\big\}\) as follows. Let \(\widehat{\mathfrak p}_{(\ell)}^{i_0}\) denote the \(\ell\)th smallest element of the sorted values of \(\big\{\widehat{\mathfrak p}_{i_0,j}: j \in {\cal I}_{\rm test}^{i_0} \big\}\). Given any \(\alpha_{\mathrm{BH}}>0\), define \[\begin{align} \widehat\ell := \max\big\{ \ell: \widehat{\mathfrak p}_{(\ell)}^{i_0} \leq \alpha_{\mathrm{BH}}\cdot \ell / |{\cal I}_{\rm test}^{i_0}| \big\}. \label{equation::BH-threshold-l} \end{align}\tag{6}\] We reject \(H_{i_0,j}\) for all \(j\) such that \(\widehat{\mathfrak p}_{i_0,j} < \widehat{\mathfrak p}_{(\widehat\ell)}^{i_0}\), and denote the rejection set by \(\mathcal{R}_{\rm BH}^{i_0}(\alpha_{\mathrm{BH}})\). The following result establishes FDR control for this local testing procedure.

2.3 Construction of \(\widehat {\mathbf{A}}\) and local tests with FDR control↩︎

To build an operational method, we construct \(\widehat {\mathbf{A}}\) that satisfies 5 for implementing the local test in Section 2.2. Due to the unknown missingness, directly applying popular existing link prediction methods, such as matrix factorization or node similarity-based tools, would fail the required exchangeabililty condition in 5 . To better illustrate our strategy, we start with an “oracle method”. Consider a weighed link estimation: \[\label{eq:weighted-est} \widehat A_{i_0,j_1} = \sum\nolimits_{j_2 \in {\cal I}_{\rm train}^{i_0}} w(j_1,j_2) A_{i_0,j_2}\tag{7}\] for \(j_1\in {\cal I}_{\rm calib}^{i_0}\cup{\cal I}_{\rm test}^{i_0}\), where the weight \(w(j_1,j_2)\) may be data-dependent. Inspired by [13], for any \(j_1 \in {\cal I}_{\rm test}^{i_0}\cup {\cal I}_{\rm calib}^{i_0}\) and \(j_2 \in {\cal I}_{\rm train}^{i_0}\), a common choice of \(w(j_1,j_2)\) is \[\begin{align} w(j_1,j_2) = \frac{K\big(d_{j_1,j_2}\big)}{\sum_{j_2' \in {\cal I}_{\rm train}^{i_0}} K\big(d_{j_1,j_2'}\big)}, \label{method::A-estimate-oracle} \end{align}\tag{8}\] where \(K(\cdot)\) is a smooth and Lipschitz kernel function, and \({d}(j_1,j_2) := \int_0^1 \Big|\int_0^1 \big\{f(u,\xi_{j_1})-f(u,\xi_{j_2})\big\} f(u,v) {\mathrm{d}}u \Big| {\mathrm{d}}v\) measures the dissimilarity of graphon slices \(f(\cdot, \xi_{j_1})\) and \(f(\cdot, \xi_{j_2})\). Since the weights \(w(j_1,j_2)\) are independent of \({\mathbf{M}}\), the unknown missing pattern does not affect the exchangeability of the entries in \(\{A_{i_0,j} - \widehat A_{i_0,j}: j \in {\cal I}_{\rm calib}^{i_0}\cup {\cal I}_{\rm test}^{i_0}\}\). As a result, the validity of FDR control for the local test via the downstream BH procedure follows naturally from Theorem 2.

To make 7 operational in practice, we replace the oracle dissimilarity \(d\) in 8 with a computable \(\widehat d\). Noted by [5], for each \(j\), we can view \(A_{i_0,j}\) as the response and \(A_{(-i_0),j}\), the remainder of column \(j\), as the predictor. The challenge, however, lies in constructing a suitable surrogate for the unobservable graphon slices \(f(\cdot, \xi_j)\)’s in 8 , such that \(\widehat d\) would not break exchangeability – for instance, simply filling the unobserved entries with 0 and plug in Equation (9) of [13] will not work.

Our key idea is to select a carefully chosen subset of \([{\mathbf{A}}\circ {\mathbf{M}}^c]_{(-i_0),j}\)’s that maintains exchangeability. This is another layer of divide-and-conquer approach. Rather than using the entire calibration set \({\cal I}_{\rm calib}^{i_0}\) to test all hypotheses in \({\cal I}_{\rm test}^{i_0}\), we divide \({\cal I}_{\rm calib}^{i_0}\) into \(|{\cal I}_{\rm test}^{i_0}|\) smaller subsets, and randomly assign one subset \({\cal I}_{\rm calib}^{i_0}(j_0)\) to each \(j_0\in{\cal I}_{\rm test}^{i_0}\). However, if \(|{\cal I}_{\rm test}^{i_0}|\) is not small relative to \(|{\cal I}_{\rm calib}^{i_0}|\), there may be insufficient calibration samples. In such nontrivial cases, we partition \({\cal I}_{\rm test}^{i_0}\) into disjoint subsets \(\{{\cal I}_{\rm test}^{i_0}(1), {\cal I}_{\rm test}^{i_0}(2), \ldots, {\cal I}_{\rm test}^{i_0}(k)\}\) and reuse the same calibration set \({\cal I}_{\rm calib}^{i_0}\) across each subset \({\cal I}_{\rm test}^{i_0}(k_0)\), for \(k_0 \in [k]\). The results are then aggregated using the e-BH procedure to form an aggregated test, as detailed in Section 2.4. For the remainder of this paper, we focus on the nontrivial case where \(|{\cal I}_{\rm test}^{i_0}|\) is not small compared to \(|{\cal I}_{\rm calib}^{i_0}|\). We therefore test hypotheses within a partition \({\cal I}_{\rm test}^{i_0}(k_0)\): \[\begin{align} H_{i_0,j}: A_{i_0,j} \leq c_{i_0,j}, ~~ \rm where j \in {\cal I}_{\rm test}^{i_0}(k_0). \label{def::hypothesis-testing-local-sampleV} \end{align}\tag{9}\] Let \(r_0\) denote the minimum size of each calibration split \({\cal I}_{\rm calib}^{i_0}(j_0)\). We choose the partition so that \(|{\cal I}_{\rm test}^{i_0}(k_0)|\leq r_1:= \lfloor |{\cal I}_{\rm calib}^{i_0}|/r_0 \rfloor\), which requires \(\lceil|{\cal I}_{\rm test}^{i_0}|/r_1\rceil\) partitions. Consequently, provided that \(|{\cal I}_{\rm calib}^{i_0}|\ge r_0\), the calibration set can be partitioned within each local test so that every test point is assigned a calibration subset of size at least \(r_0\).

With this setup, we illustrate the core of our method in Figure 2. Recall that row \(i_0\) is partitioned into training, calibration, and test sets, as shown in Figure 2-(a), and this partition aligns with the splitting strategy illustrated in Figure 1. Figure 2-(b) presents a simplified case where the calibration set is \({\cal I}_{\rm calib}^{i_0}(j_0)=\{5,6\}\), and the test set contains a single entry \(\{j_0\}:=\{7\}\). When \(|{\cal I}_{\rm calib}^{i_0}(j_0) \cup \{j_0\}|\) is small, with a high probability, we can identify multiple rows \(i\in {\cal I}_{\rm train}^{i_0}\), collectively denoted as \(\Omega_{j_0}\), for which the submatrix \({\mathbf{A}}_{\Omega_{j_0}\times ({\cal I}_{\rm calib}^{i_0}(j_0) \cup \{j_0\})}\) contains no missing entries, highlighted in dark yellow in Figure 2-(b).

Note that the construction of \(\Omega_{j_0}\) depends on the missing pattern. If \(\Omega_{j_0}=\emptyset\), we can simply define \(\widehat{d}_{j_1,j_2,j}=1\) and thus \(\widehat{d}_{j_1,j_2}=1\). It is clear that the exchangeability still holds for \(\widehat{A}_{i_0,j}\) under such uniform weight construction.

Next, we construct \(\widehat A_{i_0, j_1}\) for all \(j_1\in {\cal I}_{\rm calib}^{i_0}(j_0) \cup \{j_0\}\). Recall that we estimate \(A_{i_0,j_1}\) in 7 using a weighted average of \(A_{i_0,j_2}\)’s for \(j_2\in {\cal I}_{\rm train}^{i_0}\), where the weights depend on a dissimilarity measure \(\widehat d_{j_1,j_2}\). We now formally define \(\widehat d_{j_1,j_2}\). For each \(j\in {\cal I}_{\rm train}^{i_0}\backslash\{j_2\}\), let \(\Omega_{j_1,j_2,j}\) be the largest subset of \([n]\) satisfying: (1) \(\Omega_{j_1,j_2,j}\subseteq \Omega_{j_0}\); and (2) the submatrix \(A_{\Omega_{j_1,j_2,j}\times \{j_1,j_2,j\}}\) contains no missing entries. Let \(\langle\cdot,\cdot\rangle\) denote the inner product of two vectors. We then define \(\widehat d_{j_1,j_2,j}\) and use it to construct \(\widehat d_{j_1,j_2}\), as follows: \[\begin{align} \widehat d_{j_1,j_2,j} :=&~ \frac{\big| \langle A_{\Omega_{j_1,j_2,j}, j_1} - A_{\Omega_{j_1,j_2,j}, j_2}, A_{\Omega_{j_1,j_2,j}, j} \rangle \big|}{|\Omega_{j_1,j_2,j}|}, \tag{10} \\ \rm and\quad \widehat d_{j_1,j_2} :=&~ \sum\nolimits_{j \in {\cal I}_{\rm train}^{i_0}\backslash \{j_2\}} \frac{\widehat d_{j_1,j_2,j}}{|{\cal I}_{\rm train}^{i_0}| - 1}. \tag{11} \end{align}\] Then similar to 8 , we can construct, for \(j_1\in {\cal I}_{\rm calib}^{i_0}(j_0)\cup\{j_0\}\), \[\label{method::A-estimate} \widehat A_{i_0,j_1}:= \left\{\sum\nolimits_{j_2' \in {\cal I}_{\rm train}^{i_0}} K\big(\widehat d_{j_1,j_2'}\big)\right\}^{-1}{\sum\nolimits_{j_2' \in {\cal I}_{\rm train}^{i_0}} K\big(\widehat d_{j_1,j_2'}\big) A_{i_0,j_2'}}.\tag{12}\] For ease of reference, Table 1 summarizes the index notation introduced thus far.

Let us explain why the \(\widehat A_{i_0,j_1}\)’s constructed in such a way are exchangeable. By Proposition 1, it suffices to show that the dissimilarities \(\widehat d_{j_1,j_2}\)’s are exchangeable for all \(j_1\in {\cal I}_{\rm calib}^{i_0}(j_0) \cup \{j_0\}\), which in turn follows from the exchangeability of the underlying \(\widehat d_{j_1,j_2,j}\)’s for these \(j_1\)’s. Now we make a crucial observation that by construction, for each fixed pair \((j_2, j)\), all indices \(j_1\in {\cal I}_{\rm calib}^{i_0}(j_0) \cup \{j_0\}\) share a common \(\Omega_{j_1,j_2,j}\). For clarity, this observation is illustrated by comparing the three vertically stacked panels in Figure 2-(c), where dark orange cells highlight the rows corresponding to \(\Omega_{j_1,j_2,j}\). By Assumption 1 and model 2 , conditional on \({\mathbf{M}}\) and \(\{\xi_{j'}: j'\in {\cal I}_{\rm train}^{i_0}\}\), \(\xi_{j_1}\)’s are independent and identically distributed for \(j_1\in {\cal I}_{\rm calib}^{i_0}(j_0) \cup \{j_0\}\). This immediately implies the exchangeability of \(\widehat d_{j_1,j_2,j}\)’s, and hence that of \(\widehat A_{i_0,j_1}\)’s, for all \(j_1\in {\cal I}_{\rm calib}^{i_0}(j_0) \cup \{j_0\}\). Consequently, the nonconformity scores \(S(\widehat A_{i_0,j'}; A_{i_0,j'})\) for \(j'\in {\cal I}_{\rm calib}^{i_0}(j_0)\cup\{j_0\}\) are also exchangeable. Therefore, for \(j_0 \in {\cal I}_{\rm test}^{i_0}(k_0)\), \[\label{eqn-def::emp-p-value} \begin{align} \widehat{\mathfrak p}_{i_0,j_0} :=&~ \displaystyle\frac{ 1 + \sum_{j' \in {\cal I}_{\rm calib}^{i_0}(j_0)} \mathbb{1}\{ S(\widehat A_{i_0,j'}; A_{i_0,j'}) < S(\widehat A_{i_0,j_0}; c_{i_0,j_0})\} }{ 1 + |{\cal I}_{\rm calib}^{i_0}(j_0)| } \end{align}\tag{13}\] are valid conformal \(p\)-values [9], [21]. Moreover, they are conditionally independent by Proposition 3, making them suitable for the BH procedure.

Building on the above property of \(\widehat{\mathfrak p}_{i_0,j_0}\)’s, we apply the BH procedure in 6 to obtain a rejection set \(\mathcal{R}_{\rm BH}^{i_0,k_0}(\alpha_{\mathrm{BH}})\), as summarized in Algorithm 3. The Algorithm accommodates extreme missing patterns (e.g., no observations in an entire row or no matched indices in \(\Omega_{j_0}\)) by setting \(\widehat A_{i_0,j_1}=1\) when \({\cal I}_{\rm train}^{i_0}=\emptyset\), and by applying 12 with \(\widehat{d}_{j_1,j_2}=1\) when \(\Omega_{j_0}=\emptyset\). We conclude this section with the FDR control guarantee for Algorithm 3 in Theorem 4.

2.4 Step 2: aggregating local tests for FDR control↩︎

Recall that we break down the original problem in two stages: from an aggregated test to local tests, one for each row; then, within each row, we further break down the local test by partitioning \({\cal I}_{\rm test}^{i_0}\) into disjoint subsets. We now aggregate the results of the local tests to devise an aggregated rejection set with valid FDR control for predicting all missing links. To this end, we leverage the e-BH procedure [22], [24].

The e-BH procedure offers a notable advantage that its validity requires only that each individual e-value satisfies \(\mathbb{E}(e_{i_0,j})\leq 1\) under the null hypotheses in 9 , without any assumptions on their dependence structure. The central idea of our aggregated test is to rescale the binary test decision functions from the local tests into valid e-values. Using the leave-one-out analysis of the BH procedure [39], we show that \(\mathbb{E}\big[ \mathbb{1}\{ j\in \mathcal{R}_{\rm BH}^{i_0,k_0}(\alpha_{\mathrm{BH}}) \}/\{|\mathcal{R}_{\rm BH}^{i_0,k_0}(\alpha_{\mathrm{BH}})| \vee 1\}\big] \le \alpha_{\mathrm{BH}}/|{\cal I}_{\rm test}^{i_0}(k_0)|\), for each \(j\). Therefore, for \(j \in {\cal I}_{\rm test}^{i_0}(k_0)\), \[\begin{align} e_{i_0,j} := \frac{|{\cal I}_{\rm test}^{i_0}(k_0)| \cdot \mathbb{1}\{ j\in \mathcal{R}_{\rm BH}^{i_0,k_0}(\alpha_{\mathrm{BH}}) \}}{\{|\mathcal{R}_{\rm BH}^{i_0,k_0}(\alpha_{\mathrm{BH}})| \vee 1\}\cdot \alpha_{\mathrm{BH}}} \label{equation::compute-e-value} \end{align}\tag{14}\] is a valid e-value, namely, \(\mathbb{E}(e_{i_0,j})\leq 1\) under the null \(H_{i_0,j}\) in 4 . Then, by separately applying the e-BH procedure to \(\{e_{i_0,j}, j \in {\cal I}_{\rm test}^{i_0}(k_0)\}\), one could obtain a local test for each row \(i_0\) with the FDR controlled at level \(\alpha_{\rm BH}\). But here, we would directly construct an aggregated test using all \(\{e_{i,j}:(i,j) \in {\cal D}_{\rm test}\}\). To this end, we sort all e-values \(\{e_{i,j}:(i,j) \in {\cal D}_{\rm test}\}\) in descending order: \(e_{(1)} \geq \cdots \geq e_{(|{\cal D}_{\rm test}|)}\). Then, for any given \(\alpha_{\mathrm{eBH}}>0\), set \[\begin{align} \widehat k := \max\left\{ k: 1\leq k\leq |{\cal D}_{\rm test}|, e_{(k)} \geq \frac{|{\cal D}_{\rm test}|}{\alpha_{\mathrm{eBH}}k} \right\}~ \rm or ~ \widehat k = 0 \rm if the set is empty, \label{eqation::compute-hat-k} \end{align}\tag{15}\] and reject all \(H_{i,j}\)’s with e-values \(e_{i,j}\geq e_{(\widehat k)}\). The procedure fails to reject any hypotheses when \(\widehat k = 0\). This rejection rule controls FDR of 1 at the level \(\alpha_{\mathrm{eBH}}\).

To mitigate variations arising from the random selection of \({\cal I}_{\rm calib}^{i_0}\) and \({\cal I}_{\rm train}^{i_0}\), we further adopt a derandomization procedure [24]. This builds on the property that the average of multiple e-values remains a valid e-value [23]. For each row \(i_0\) and \(j \in {\cal I}_{\rm test}^{i_0}(k_0)\), we repeat Algorithm 3 for \(m_{i_0}\) times and compute the e-values \(e_{i_0,j}^{(1)},\ldots, e_{i_0,j}^{(m_{i_0})}\). Then we derandomize the e-values as follows: \[\begin{align} \bar{e}_{i_0,j} := \frac{1}{m_{i_0}} \sum_{l=1}^{m_{i_0}} e_{i_0,j}^{(l)}, ~~ \rm wherej \in {\cal I}_{\rm test}^{i_0}(k_0). \label{equation::compute-e-value-avg} \end{align}\tag{16}\] Finally, we apply the e-BH procedure to all \(\bar e_{i,j}\)’s as in 15 , using them in the same manner as \(e_{i,j}\) in the original algorithm described above. This yields our method clp, summarized in Algorithm 4, which enjoys finite-sample valid FDR control as shown in Theorem 5.

We make two key remarks. First, our method clp accommodates arbitrary missing mechanism, provided Assumption 1 hold. Examples include: (1) random but non-uniform missingness [27], and (2) structured missingness, such as block or staggered patterns commonly encountered in panel data [40], [41]. Second, as shown in Theorem 5, clp guarantees valid FDR control for any choice of \(\alpha_{\mathrm{BH}}\), while [24] recommend setting \(\alpha_{\mathrm{BH}}=\alpha_{\mathrm{eBH}}/2\) for improved powers.

We note that our method may exhibit some conservativeness in practice. By 14 , the maximum possible value of \(e_{i_0,j}\) for \(j \in {\cal I}_{\rm test}^{i_0}(k_0)\) is \(|{\cal I}_{\rm test}^{i_0}(k_0)|/\alpha_{\mathrm{BH}}\). Thus, when \(|{\cal I}_{\rm test}^{i_0}(k_0)|\) is small due to extensive data-splitting, the resulting e-values may be lowered, limiting the number of rejections under the e-BH procedure in 15 . To alleviate this effect, we introduce an inflation factor to adjust each e-value \(e_{i,j}\) and increase the likelihood of rejection. Numerical studies in Section 3 explore distinct choices of this inflation factor and suggest that using a factor of \(1/\alpha_{\mathrm{BH}}\) effectively enhances power while maintaining well-controlled FDR. A more rigorous theoretical characterization of this practically effective factor remains an intriguing yet challenging direction for future work.

Computing the dissimilarity measure in 11 is computationally intensive due to the evaluation of all terms \(\widehat d_{j_1,j_2,j}\) in 10 . The overall computational complexity of clp is \(O\big(m|{\cal I}_{\rm train}^{i_0}|^2|{\cal I}_{\rm calib}^{i_0}|\big)\), where \(m\) is the number of missing entries, assuming \(|{\cal I}_{\rm train}^{i_0}|\) (and \(|{\cal I}_{\rm calib}^{i_0}|\)) are of the same order across all \(i_0 \in [n]\). Specifically, computing 10 and 11 incurs a cost of \(O\big(|{\cal I}_{\rm train}^{i_0}|\big)\), while the aggregation step in 12 adds additional \(O\big(|{\cal I}_{\rm calib}^{i_0}|)\). In practice, our method can be accelerated by reducing the sizes of \({\cal I}_{\rm train}^{i_0}\) and \({\cal I}_{\rm calib}^{i_0}\); for instance, setting \({\cal I}_{\rm train}^{i_0}\asymp{\cal I}_{\rm calib}^{i_0}=O(1)\) reduces the overall complexity to \(O(m)\), though at the potential cost of reduced accuracy in estimating dissimilarity and, consequently, lower power.

2.5 Extensions and comparisons↩︎

While our proposed framework is detailed for directed networks, it naturally extends to bipartite and undirected networks with minimal modifications. The adaptation to bipartite networks is straightforward, owing to the their exchangeability. That is, we can use more rows to construct \(\widehat A_{i_0,j}\) and the associated \(p\)-values. Extending to undirected networks requires extra care due to the symmetry constraints \(A_{i,j}=A_{j,i}\) and \(M_{i,j}=M_{j,i}\), but the solution is remarkably simple. We first apply our Algorithms 3 and 4 up to Line 9 of Algorithm 4, where the \(\bar e_{i,j}\) values are obtained. Thereafter, we apply the remainder of Algorithm 4, the e-BH procedure, only to the upper triangular part, i.e., \(\{\bar e_{i,j}: (i,j)\in {\cal D}_{\rm test}\rm and i<j\}\), instead of the full set \(\{\bar e_{i,j}: (i,j)\in {\cal D}_{\rm test}\}\). Formal algorithmic descriptions (Algorithms 16 and 18) and theoretical justifications for both extensions are provided in Section 8 of the Supplementary Material.

We conclude this section by comparing our method with existing works. [27] considers arbitrary \({\mathbf{A}}\) under a missing-at-random mechanism but requires knowledge of, or accurate estimates for, the odds ratios introduced by missing probabilities. In contrast, we allow arbitrary missingness patterns \({\mathbf{M}}\) that are independent of \({\mathbf{A}}\), and leverage the latent graphon structure in \({\mathbf{A}}\) to reinstate the exchangeability. This enables our method to accommodate a broader class of missing mechanisms. Overall, the two approaches rely on different assumptions and do not subsume each other. We also note that our paper considers random matrices \({\mathbf{A}}\), whereas in related work such as [27], [30], the matrix to be tested is treated as fixed. More recently, [28], [29] proposed a novel FDR control approach that casts link prediction as a classification problem. While innovative, their formulation is inherently limited to unweighted networks and assumes missingness rates take only two distinct values, depending on whether an entry is \(0\) or \(1\). In contrast, our method applies to both weighted and unweighted networks and allows fully heterogeneous, unknown missingness rates, offering greater flexibility in practice. It is also possible to leverage general p-to-e calibrators [23], [42] to construct the e-values, but as reported in literature, this approach tends to be rather conservative. We leave further studies along this line to future research.

3 Simulations↩︎

3.1 Simulation 1: Tuning our algorithm↩︎

We generate networks following the distribution \(A_{i,j} = f(\xi_i,\xi_j) + \epsilon_{i,j}\), where \(\epsilon_{i,j} \sim \mathrm{Uniform}[-0.1,0.1]\). Three graphons are considered to cover a wide range of scenarios:

(i) Setting 1 (low-rank and smooth): \(f_1(\xi_i,\xi_j) = \xi_i^3 + 2\xi_j^3\);

(ii) Setting 2 (high-rank): \(f_2(\xi_i,\xi_j) = \max(\xi_i, \xi_j)^{\frac{2}{3}}\cos[0.1/\{(2 \xi_i-1/2)^3 + (\xi_j-1/2)^3 + 0.01\}]\);

(iii) Setting 3 (non-smooth): \(f_3(\xi_i,\xi_j) = (3 \xi_i^2 + \xi_j^2) \cdot \cos\{{(2 \xi_i^4 + \xi_j^4)}^{-1}\}\).

We fix the network size at \(n=200\) and set the proportion of missing links to \(0.2\), resulting in an average of \(3,980\) missing links. Specifically, we generate \(M_{i,j}\sim\) Bernoulli\((q_{i,j})\) and consider heterogeneous missing probabilities, where \(q_{i,j}\)’s are independently generated from \(\mathrm{Uniform}[0,0.4]\). We set the null hypotheses to be true, namely \(c_{i,j} = A_{i,j}\) in 1 , for \(70\)% of entries. For the remaining \(30\%\), we set \(c_{i,j} = A_{i,j} - 1.5\). For each row \(i_0\), all observed entries are randomly split into \({\cal I}_{\rm train}^{i_0}\cup {\cal I}_{\rm calib}^{i_0}\) with a ratio of \(|{\cal I}_{\rm train}^{i_0}|:|{\cal I}_{\rm calib}^{i_0}|=2:3\), and we set \(\alpha_{\mathrm{BH}}= \alpha_{\mathrm{eBH}}/2\) following [24]. In the BH procedure, we choose \(r_0=25\), as defined in the paragraph following 9 . We use the Gaussian kernel \(K(x) = \exp\big(-x^2/2\big)/\sqrt{2\pi}\) in 12 and derandomize local tests using \(20\) repetitions in 16 .

Recall that in 14 , we inflate the e-values by a factor of \(c/\alpha_{\mathrm{BH}}\). Here, we examine the impact of \(c\in\{0.5, 1, 1.5\}\). Performance is evaluated using the FDR, approximated by the averaged FDP, and power, both assessed over \(100\) experiments. We vary \(\alpha_{\mathrm{eBH}}\) from \(0.1\) to \(0.3\). The results in Figure 5 show that our method consistently controls the empirical FDR across significance levels when \(c=0.5\), and in most settings when \(c=1\). It also achieves desirable power results as \(\alpha_{\mathrm{eBH}}\) increases. Overall, we recommend \(c=1\) as it strike a satisfactory balance between FDR control and decent power.

3.2 Simulation 2: Comparison with benchmark methods↩︎

We compare our method (clp) with conformal matrix completion (cmc, [27]) and the link prediction method by [28] (dss) via two sub-simulations for weighted and unweighted networks, respectively. While cmc applies to both settings, dss is, by design, primarily suited for unweighted networks.

For the first sub-simulation (weighted links), we adopt Settings 1–3 from Section 3.1, fixing \(\alpha_{\mathrm{BH}}=0.1\) and \(\alpha_{\mathrm{eBH}}=0.2\). The threshold \(c_{ij}\) is set as a specific quantile of the observed link weights, with the quantile level \(\kappa\) varying from \(0.1\) to \(0.5\). As noted earlier, we inflate the e-values in 14 by a factor of \(1/\alpha_{\rm BH}\). Note that cmc only provides confidence intervals with guaranteed coverage, rather than \(p\)-values. Its underlying idea is to approximate the population quantiles of absolute prediction errors by the empirical quantiles of calibrated absolute errors, a procedure valid only for two-sided problems. To adapt it to our one-sided multiple testing problem 1 , we modify it to produce comparable \(p\)-values: we omit the absolute value when computing errors and construct the empirical distribution function \(\widehat{F}_{\rm eCDF}(\cdot)\) based on the calibrated errors \(\widehat{A}_{ij}-A_{ij}\). We then compute \(p\)-values for cmc as \(1-\widehat{F}_{\rm eCDF}(\widehat{A}_{ij}-c_{ij}), (i,j)\in{\cal D}_{\rm test}\), and apply the BH procedure. Note also that cmc requires a hypothesized rank as input. In practice, results of cmc are robust to the hypothesized rank, as also noted by [27]. Accordingly, we fix the rank at \(r=2\) for cmc across all settings for demonstration.

As shown in Figure 6, while both clp and cmc control the FDR, clp consistently achieves higher power across all settings. We note that cmc is sensitive to the choice of \(c_{i,j}\) values. In Section 6 of the Supplementary Material, we adopt the approach in Section 3.1 to select \(c_{i,j}\) based on signal strength. Under such configuration, cmc fails to control the FDR, whereas clp continues to perform well.

To compare with dss, we set up a second sub-simulation for unweighted networks, using a different network-generating process from that used for the weighted case. Specifically, we generate \(A_{i,j} = \mathbb{1}\{(\xi_i+\xi_j)/2 + \epsilon_{i,j}>t\}\), where \(\epsilon_{i,j} \sim \mathrm{Uniform}[-0.25,0.25]\). We set \(c_{i,j}=0.5\) in 1 and vary \(t \in \{0,0.05,0.1,0.15,0.2\}\) to control the number of null and alternative hypotheses. We consider heterogeneous missingness by drawing \(q_{i,j}\sim \mathrm{Uniform}[0,0.2]\) independently for each \((i,j)\), yielding an average proportion of missing links of approximately \(10\%\). We fix \(\alpha_{\mathrm{BH}}=0.2\) and \(\alpha_{\mathrm{eBH}}=0.2\), and use \(r_0=50\) to improve estimation of \(A_{i,j}\) for unweighted networks.

As displayed in Figure 7, all three methods control the FDR across different values of \(t\). The conservativeness of cmc aligns with our observation for weighted networks. While the FDR of dss exceeds that of clp, its relative low power can be attributed to the heterogeneous missing probabilities. In general, both clp and dss outperform cmc in the unweighted network setting, and our method demonstrates strong performance across both weighted and unweighted networks.

4 Data example: bilateral trade flows↩︎

We consider the bilateral trade flow data from [10], which records trade volumes between \(158\) countries in 1989. The dataset forms a directed network, where both indices \(i\) and \(j\) represent countries acting as exporters and importers, respectively. The outcome \({A_{i,j}}\) denotes the log-transformed trade volume from country \(i\) to country \(j\), with a mean of \(3.82\), standard deviation of \(4.69\), and a range from \(0\) to \(18.39\). Here, \(A_{i,j}=0\) indicates no recorded trade flow. We randomly mask \(10\%\) of the links, corresponding \(2,496\) trade entries, for evaluating the prediction performance of different methods. For these held-out entries, we set \(c_{i,j} = 0.2\) in the random hypotheses 1 , resulting in approximately equal proportion of null and alternative hypotheses.

We compare our method only with cmc, as the method dss is not applicable because it is classification-based and limited to binary or categorical links. Following our simulation setup, we implement clp with \(|{\cal I}_{\rm train}^{i_0}|:|{\cal I}_{\rm calib}^{i_0}|=2:3\) for each \(i_0 \in [n]\), set \(r_0=25\), and use \(\alpha_{\mathrm{BH}}=\alpha_{\mathrm{eBH}}/2\). We employ a Gaussian kernel for \(K(\cdot)\) in 8 and perform \(20\) derandomizations in 16 . The nominal FDR level \(\alpha\) is varied over \(\{0.1, 0.2, 0.3\}\).

The left panel of Figure 8 shows that clp controls FDR well across all the nominal FDR levels, whereas cmc substantially inflates the FDP. The right panel demonstrates that our method is also more powerful than cmc. As \(\alpha\) increases, the empirical power of clp improves, and its advantage over cmc becomes more pronounced.

5 Concluding remarks↩︎

In this article, we develop a general framework for FDR control in multiple network link prediction, leveraging e-values and insights from conformal inference. We propose data-driven methods and a multiple-splitting strategy to reinstate the exchangeability required for computing nonconformity scores under complex network edge dependencies and heterogeneous missingness. With supporting theories, we establish performance guarantees for multiple link prediction without requiring knowledge of the underlying network distribution.

Experiments on both simulated and real data suggest that a \(1/\alpha_{\mathrm{BH}}\) inflation factor performs well for clp. Setting a conservative threshold for e-values in the e-BH procedure may be inefficient in practice. [43] proposed some preliminary strategies for boosting e-values in the e-BH procedure under certain distributional assumptions while maintaining FDR control. It is of interest to investigate whether an optimal inflation criterion of e-values exists in our context, and if so, how it can be effectively designed.

We note that block-wise schemes may be designed to offer computational savings. However, such approaches are challenging in our setting, as the required block-level exchangeability is fragile under random missingness in \({\mathbf{M}}\). Missing entries can break the symmetry needed for valid permutation or aggregation, making it difficult to construct blocks that both satisfy exchangeability and remain robust to heterogeneous missingness patterns. Moreover, forming blocks using only observed entries limits their applicability for testing missingness itself. We leave a rigorous treatment of block-wise aggregation under missingness to future work.

Code↩︎

The code for reproducibility is available at: https://github.com/code-cloud9/CLP.

Acknowledgement↩︎

Dong Xia was supported by Hong Kong RGC grant GRF 16303224. Yuan Zhang was supported by United States NSF grant DMS-2311109. Wen Zhou was supported by United States NIH grants R01GM157600 and R01GM144961.

This Supplementary Material contains proofs of Theorems 2–5 and Propositions 1–3 in Section A. Additional results of numerical studies are provided in Section 6 for simulations and Section 7 for real data application. The extensions of our proposed method clp to undirected and bipartite networks are given in Section 8. Additional results for Bernoulli missing patterns are given in Section 9. Unless stated otherwise, all the notation follows the same definitions as in the main body of the paper.

6 Additional numerical studies↩︎

We explore how signal strength affect the performance of clp, cmc and dss. For weighted network setting, we apply the graphons in Section 3.1, setting \(c_{i,j} = A_{i,j} - \delta\) for \(30\)% of the hypotheses and set \(c_{i,j} = A_{i,j}\) for the rest as the null hypothesis. Here, \(\delta\) controls signal strength.

Figure 9 shows the empirical FDR and power across different values of signal strength \(\delta\). As can be seen, cmc fails to control the FDR across all settings. Comparing the results in Figure 9 and Figure 6, we see that the performance of cmc is sensitive to the choice of \(c_{ij}\)’s. In contrast, clp consistently achieves desirable FDR control and its power increases as the signal strength becomes larger across all settings.

Since the signal strength is not straightforward to define in the unweighted setting, we slightly modify our testing target to illustrate how signal strength affects the performance of different methods. Specifically, we generate \(p^*_{i,j} = f(\xi_i,\xi_j) + \epsilon_{i,j}\) with the graphon functions in Settings 1–3 and rescale them by \(p_{i,j}=(p^*_{i,j}-\min_{i,j}p^*_{i,j})/(\max_{i,j}p^*_{i,j}-\min_{i,j}p^*_{i,j})\) to \([0,1]\)-valued probabilities. Then we independently generate \(A_{i,j} \sim \rm Bernoulli(p_{i,j})\). Different from the previous settings, here we consider heterogeneous missing probabilities where \(q_{i,j}\)’s are independently generated from \(\mathrm{Uniform}[0,0.4]\). Instead of testing whether there is a true link or not, we test \(H_{i,j}: p_{i,j} \leq c_{i,j}\) for the random \(p_{i,j}\)’s, and set \(c_{i,j} = p_{i,j} - \delta\) for 30% of the hypotheses regarding the missing entries. Otherwise, we set \(c_{i,j} = p_{i,j}\) for the null hypotheses.

As shown in Figure 10, clp and dss control FDR throughout the three settings, while cmc fails. Additionally, the relative low power of dss compared to clp can be attributed to the heterogeneous missing probabilities. [28] assumes that the missingness rates \(q_{i,j}\)’s take only two distinct values, depending on whether the links are \(0\) or \(1\). The target of testing random \(p_{i,j}\)’s instead of \(A_{i,j}\)’s can also contribute to dss making fewer rejections.

6.1 Comparison with dss under weighted networks↩︎

Note that the procedure in [28] (dss) cannot be easily extended to weighted network settings with heterogeneous null values \(c_{i,j}\), since their \(D_{\rm train}^0\) (the set of observed non-existent edges) and \(D_{\rm train}^1\) (the set of observed true edges) are not well-defined when \(c_{i,j}\) differs across missing entries. Furthermore, in Settings \(1\)–\(3\) introduced in Section 3.1 of the main paper, where the \(c_{i,j}\)’s are uniform, one way to apply the procedure in [28] to test \(H_{i,j}: A_{i,j} \leq c_{i,j}~ \rm for (i,j) \in {\cal D}_{\rm test},\) is as follows: during testing, we first transform the observed weighted network into an unweighted network using the threshold \(c_{i,j}\) (set to the quantile of the observed links in Settings \(1\)–\(3\)). We then apply the procedure in [28] to this transformed unweighted network.

The results are shown in Figure 11. We observe that the procedure in [28] fails to control the FDR. This can be attributed to the heterogeneous missing probabilities in Settings \(1\)–\(3\), which are not addressed by their method, as it assumes homogeneous missing probabilities.

6.2 Sensitivity of cmc to the choice of hypothesized rank↩︎

For the cmc method, we conducted additional experiments using a larger hypothesized rank, namely \(r=10\), under Settings 1–3 introduced in Section 3.1 of the main paper. As shown in Figure 12, the performance of cmc with \(r=10\) is similar to that with \(r=2\) in these settings, in terms of both empirical FDR and power. As noted by [27], in practice, results of cmc are robust to the hypothesized rank. Therefore, the performance of cmc shown in Figure 12 can be attributed to the inherent robustness of the cmc method and that in Settings \(1\)–\(3\), increasing \(r\) to \(10\) may not alter the ability to recover the core low-dimensional signal in the network.

6.3 Weighted and unweighted estimators in clp↩︎

Recall that in Section 2.3 of the main paper, we construct a “fine-grained weighted estimator" of \({\mathbf{A}}\) through the weights \(w(j_1,j_2),\) which are determined by a dissimilarity measure \(d(j_1,j_2)\). The purpose of this construction is to preserve the relevant exchangeability structure while improving the estimation accuracy of \(\widehat{{\mathbf{A}}}\). A natural question is that whether one could instead use an unweighted estimator, that is, an estimator with constant weights, and thereby pursue a trade-off between estimator simplicity and the estimation accuracy. The answer is yes. In fact, the unweighted estimator can be viewed as a special case of our general construction. Specifically, if we take \(d(j_1,j_2)\equiv 1,\) then the resulting weight becomes \(w(j_1,j_2)=\frac{1}{|{\cal I}_{\rm train}^{i_0}|},\) which is constant over \(j_2\in {\cal I}_{\rm train}^{i_0}\). Hence, the corresponding estimator is unweighted. Moreover, the exchangeability argument underlying our clp procedure continues to apply in this case. Therefore, from a validity perspective, the unweighted estimator is also admissible.

However, using constant weights may substantially degrade the estimation accuracy of \(\widehat{\mathbf{A}}\), which may consequently reduce power. To illustrate this trade-off, we conduct additional simulations comparing the weighted estimator and the unweighted estimator. We consider both a weighted-network setting and an unweighted-network setting. For the weighted-network setting, we use Setting 3 in Section 3.1 of the main paper. For the unweighted-network setting, we use the data-generating mechanism introduced in Section 3.2. In both cases, we apply clp with the same implementation as in the main paper, including the \(1/\alpha_{\rm BH}\) inflation factor for the e-values. Figure 13 reports the empirical FDR and power for the two estimators. We present our original clp results using the “fine-grained weighted estimator”, shown in blue and labeled as the “weighted estimator”. For comparison, the results using the “unweighted estimator" sets \(d\left(j_1, j_2\right) \equiv 1\), as discussed above. Its performance is shown in orange. Under the weighted-network setting, the weighted estimator and the unweighted estimator perform similarly. Under the unweighted-network setting, the weighted estimator yields better power than the unweighted estimator, while still maintaining FDR control. These results suggest that although the unweighted estimator is valid and conceptually simpler, the fine-grained weighted estimator can have better power performance in practice.

6.4 Sensitivity of clp to the choice of \(r_0\)↩︎

In our construction, the e-value for \(j \in \mathcal{I}_{\rm test}^{i_0}(k_0)\) contains the factor \(|\mathcal{I}_{\rm test}^{i_0}(k_0)|\) in the numerator, and the BH rejection set used in its denominator is also computed within this local test set. Therefore, the choice of \(r_1\), which controls the upper bound on \(|\mathcal{I}_{\rm test}^{i_0}(k_0)|\), can affect the magnitude of the resulting e-values and hence the empirical performance of the final e-BH procedure. Specifically, a larger size of the local test set (larger \(r_1\)) may inflate e-values and potentially increase power, but it also reduces the size of each calibration subset \(\mathcal{I}_{\rm calib}^{i_0}(j_0)\), which may lead to more conservative conformal p-values in the BH procedure. Recall that in our clp implementation, \(r_1\) is determined by the calibration split size \(|{\cal I}_{\rm calib}^{i_0}|\) and \(r_0\) through \(r_1 := \lfloor |\mathcal{I}_{\rm calib}^{i_0}| / r_0 \rfloor\). Therefore, \(r_1\) is not an entirely independent tuning parameter, but can be tuned indirectly through \(r_0\), where a larger value of \(r_0\) corresponds to a smaller \(r_1\). To examine the sensitivity of clp to the choice of \(r_0\) (and hence \(r_1\)), we conducted additional experiments under Settings 1–3 from Section 3.1 of the main paper. We set \(\kappa = 0.3\) and vary \(r_0 \in \{10, 20, 30, 40\}\). The results are shown in Figure 14.

Across all three settings, the empirical FDR remains below the nominal level for all values of \(r_0\) considered, confirming that FDR control is robust to this choice. Additionally, the empirical power in Settings 1 and 2 increases slightly as \(r_0\) increases, while in Setting 3 it remains nearly unchanged. Overall, these results suggest that the performance of our method is not sensitive to the local test set size over the considered range, and that \(r_0 = 25\) as used in the main paper provides a reasonable default.

7 Additional real data analysis↩︎

In Section 4, we set \(c_{i,j}=0.2\) and vary \(\alpha \in \{0.1, 0.2, 0.3\}\). In this section, we fix \(\alpha=0.2\) and vary the value of \(c_{i,j}\) to further demonstrate the performance of clp and cmc. Specifically, we consider \(c_{i,j} \in \{0.2, 2, 4.3, 7.1\}\), with the proportion of alternative hypotheses gradually decreasing to approximately 30%. The empirical FDR and power (averaging across \(100\) independent runs) are summarised in Figure 15. As expected, clp achieves valid FDR control across all configurations, while cmc does not. Moreover, as \(c_{i,j}\) becomes larger, the power of both methods decreases. As a result, clp performs well in detecting whether trading activity exists between two countries but becomes less powerful when detecting higher volumes of bilateral trade.

8 FDR control for undirected and bipartite networks↩︎

In this section, we extend the proposed clp framework to accommodate undirected and bipartite networks, elaborating on the modified version of Algorithm 3. Notice that once Algorithm 3 and its properties are still valid to control the FDR for the local test as defined in Theorem 4, we can directly incorporate it into Algorithm 4 and establish the FDR control across all missing links by applying the proof of Theorem 5. We start with the undirected networks. Due to symmetry, we consider that information on the network and missing patterns is incorporated in the upper triangular part, i.e., \(\{A_{i,j}\}\), \(\{M_{i,j}\}\), for \(i<j\). To facilitate the analysis, we write the extended version of \({\mathbf{A}}\), \({\mathbf{M}}\) as \({\mathbf{A}}'\), \({\mathbf{M}}'\), where \[\label{eq:A-M-extend} A_{i,j}'=\left\{\begin{array}{cc} A_{i,j} & \rm if{ i<j} \\ A_{j,i} & \rm if{ i>j} \end{array}\right. ~~; \quad M_{i,j}'=\left\{\begin{array}{cc} M_{i,j} & \rm if{ i<j} \\ M_{j,i} & \rm if{ i>j} \end{array}\right. .\tag{17}\] We now describe the FDR control for the local test in undirected networks in Alogrihtm 16.

Figure 17 provides an illustrative example of the splitting procedure in clp for an undirected network. It is seen that the construction follows the same procedure as described in Section 2.2 for directed networks. Based on Proposition 6, the FDR control results below follows directly by the proof of Theorems 4–5.

We then turn to the bipartite networks, where \({\mathbf{A}}, {\mathbf{M}}\in \mathbb{R}^{n \times m}\) and the bipartite network \({\mathbf{A}}\) follows the decomposition: \[\begin{align} A_{i,j}\sim f(\xi_i,\zeta_j), \rm ~for~i \in [n], \;j\in[m], \label{def::exch-model-bipartite} \end{align}\tag{18}\] where \(\{\xi_i\}_{1\leq i\leq n}, \{\zeta_j\}_{1\leq j\leq m}\stackrel{\rm i.i.d.} \sim \mathrm{Uniform}[0,1]\). Notice that, the bipartite networks can be viewed as a \(n\times m\) sub-network of the directed network 2 in \(\mathbb{R}^{ (n+m) \times (n+m)}\). Suppose that \({\mathbf{A}}^{\rm full}\in\mathbb{R}^{ (n+m) \times (n+m)}\) is a directed network generated by \(\big\{\{\xi_i\}_{1\leq i\leq n}, \{\zeta_j\}_{1\leq j\leq m}\big\}\stackrel{\rm i.i.d.} \sim \mathrm{Uniform}[0,1]\) following 2 . Then, it is seen that \({\mathbf{A}}={\mathbf{A}}^{\rm full}[1:n,n+1:n+m]\) is a submatrix of the exchangeable \({\mathbf{A}}^{\rm full}\). Following this observation, we can show the exchangeablity of bipartite networks below.

The poof of Lemma 1 follows the property 3 by noticing that any \(\pi_1\) on \([n]\) and \(\pi_2\) on \([m]\) can be combined as a unique \(\pi_{\rm full}\) acting on \([n+m]\) with \[\pi_{\rm full}(i)= \left\{ \begin{array}{cl} \pi_1(i) & \rm wheni\in [n], \\ \pi_2(i-n)+n & \rm wheni \in[n+m] \backslash [n]. \\ \end{array}\right.\] We now modify our Algorithm 3 to the following version by substituting \({\cal I}_{\rm train}^{i_0}\) in step 3 by \([n]{\backslash \{i_0\}}\):

By the construction of \(\Omega_{j_0}\), the submatrix \([{\mathbf{A}}\circ {\mathbf{M}}^c]_{\Omega_{j_0},{\cal I}_{\rm calib}^{i_0}(j_0)\cup\{j_0\} }\) is always fully observed. According to Lemma 1, we can show the following exchangeability property within the submatrix.

An example illustrating the splitting procedure in clp for bipartite networks is shown in Figure 19. Now we are ready to construct the conformal p-values. We can analogously define \(\Omega_{j_1,j_2,j}\) as in the main paper. By construction, for every given \((j_2, j)\), where \(j_2 \in {\cal I}_{\rm train}^{i_0}\), all \(j_1\in {\cal I}_{\rm calib}^{i_0}(j_0) \cup \{j_0\}\) share a common \(\Omega_{j_1,j_2,j}\). By Proposition 8, conditional on the training data \(\{\xi_i\}_{i\in\Omega_{j_0}}\), \(\{\zeta_j\}_{j\in{\cal I}_{\rm train}^{i_0}}\), training block \(\{A_{i,j}\}_{i\in\Omega_{j_0},j\in{\cal I}_{\rm train}^{i_0}}\) and missing pattern \({\mathbf{M}}\), it is seen that for each \(j_2\in {\cal I}_{\rm train}^{i_0}\), we have \[\begin{align} \widehat{d}(\pi(j_1),j_2) =&~ \sum_{j \in {\cal I}_{\rm train}^{i_0}\backslash \{j_2\}} \frac{ |\langle A_{\Omega_{\pi(j_1),j_2,j},\pi(j_1)} - A_{\Omega_{\pi(j_1),j_2,j},j_2}, A_{\Omega_{\pi(j_1),j_2,j},j} \rangle| }{ \big( |{\cal I}_{\rm train}^{i_0}| - 1 \big) |\Omega_{\pi(j_1),j_2,j}| } \\ \stackrel{d}{=} & \sum_{j \in {\cal I}_{\rm train}^{i_0}\backslash \{j_2\}} \frac{ |\langle A_{\Omega_{j_1,j_2,j},j_1} - A_{\Omega_{j_1,j_2,j},j_2}, A_{\Omega_{j_1,j_2,j},j} \rangle| }{ \big( |{\cal I}_{\rm train}^{i_0}| - 1 \big) |\Omega_{j_1,j_2,j}| }\\ =&~~ \widehat{d}(j_1,j_2). \end{align}\] The second equality above is due to the the definition of inner product and Proposition 8. This directly gives the following proposition, which can be justified by applying the proof of Proposition 3.

Analogous to Theorem 7, according to Proposition 9, we obtain the following FDR control results.

9 Additional theoretical results↩︎

We now show the existences of \(\Omega_{j_0}\) for matching that leads to non-trivial rejections of our algorithm. The proposition below guarantees that, when \(M_{i,j}\) are independent Bernoulli\((q_{i,j})\) random variables with unknown and potentially heterogeneous missing probabilities \(q_{i,j}\), there exists a sufficiently large number of fully observed rows to construct \({\mathbf{A}}_{\Omega_{j_0}\times ({\cal I}_{\rm calib}^{i_0}(j_0) \cup \{j_0\})}\), provided the missing rates satisfy the mild condition \(q_{i,j} \le c_0 < 1\).

Moreover, we show that under structured missingness patterns, such as block or staggered missingness [40], [41], there remain sufficiently many fully observed rows for each \((i_0,j_0)\) under mild regularity conditions. We take block missingness as an illustrative example; the argument for staggered missingness follows analogously.

References↩︎

---

1. Equal contribution.↩︎

2. Corresponding author: mailto:yzhanghf@stat.osu.edu↩︎