Flexible Imputation of Incomplete Network Data

Flexible Imputation of Incomplete Network Data

Keywords: Sampled network, network formation model, imputation, pseudo-distance, GMM

\thispagestyle{empty}

\setcounter{page}{1}

Introduction

There is a large and growing literature studying how social networks shape individuals' behavior and economic outcomes. Examples include how social interactions affect information diffusion banerjee2013diffusion, beaman2018diffusion, test scores sacerdote2001peer, and demand for financial assets conley2010learning, cai2015social. While network structure is central to these analyses, collecting complete network data is often expensive breza2020using, so applied work typically relies on sampled or partially observed networks instead. However, it has been noted that using sampled networks in place of the full network may lead to biased regression coefficients and generalized method of moments (GMM) estimators; see, for example, chandrasekhar2011econometrics. As a result, many empirical conclusions may be undermined by incomplete network data.

To overcome this limitation, we propose a flexible imputation method for sampled networks and show that downstream empirical analysis using the imputed network delivers consistent parameter estimators. Our method is designed for networks generated by dyadic formation processes and observed through egocentric sampling, a widely used design in which researchers randomly sample\footnote{ Our method also applies to more general sampling schemes in which individuals are randomly sampled with heterogeneous probabilities that depend on both observed and unobserved characteristics. } a subset of individuals and record all of their links in the network, not only those within the sample. The proposed procedure imputes missing links by combining two parts: (i) one explained by observable covariates, recovered via projection onto the observed covariates, and (ii) a residual part estimated using a local two-way fixed-effects regression. The method has several advantages over existing approaches, making it particularly well suited to empirical studies of social networks in economics.

First, the proposed method avoids parametric assumptions. Many imputation approaches, such as chandrasekhar2011econometrics, breza2020using, and boucher2020estimating, require specifying the functional form of the network formation model. In practice, this is a demanding requirement, as applied researchers rarely have reliable prior information about the correct functional form. Our approach circumvents this issue by using a fully nonparametric procedure and therefore does not require specifying the network formation process.

Second, our imputation method does not rely on a low-rank assumption. Assuming a low-rank network formation model and exploiting this structure for imputation is a common approach in the statistical network literature li2023link, cai2016structured, but this assumption may be violated when the network formation process is nonlinear, as in the specification in graham2017econometric. Building on ideas from freeman2023linear and feng2023optimal, our method instead approximates network formation using local two-way fixed-effects regressions, making it applicable to a broad class of nonlinear dyadic formation processes.

Third, our procedure jointly exploits observed covariates and accommodates unobserved heterogeneity in link formation. Social networks typically exhibit homophily---meaning that individuals with similar characteristics are more likely to form links---and degree heterogeneity---individuals differ systematically in the number of links. An effective imputation method should therefore exploit information from observed covariates and accommodate latent features. The method proposed in this paper achieves both objectives without imposing additional structure on the joint distribution of latent factors and observed covariates.

The general nonparametric network formation setting leads to a technical challenge: because the underlying formation model is unknown and potentially nonlinear, we cannot directly estimate latent factors, as is common in linear factor models, and then use them to recover the formation process and impute missing links. To address this problem, we draw on the insight of zhang2017estimating that similarity in observed connections can reveal similarity in latent factors, and use the sampled subnetwork to construct a pseudo-distance that serves as a proxy for latent distance. For each unsampled individual, we then use this pseudo-distance to locate sampled individuals who are “close" to her and use their observed connections to impute her missing links.

We establish a convergence rate for the imputed link probabilities as the network size and the number of sampled individuals go to infinity. The convergence is uniform over all missing links, which implies that the imputation error is controlled at the link level and is particularly useful for applications that require accurate imputation of each link probability. We further provide a bias-variance decomposition of the imputation error, which, to the best of our knowledge, is new in this context. The decomposition shows that the approximation error introduced by using pseudo-distance has an asymmetric effect: it is asymptotically negligible for the bias, but contributes to the variance.

We then establish the consistency of the GMM estimator constructed from the imputed network. This result is particularly relevant for empirical applications, where interest typically centers on the consistency of estimators based on the imputed network rather than on the accuracy of the imputed links themselves. In addition, because our imputation procedure does not require specifying a parametric network formation model and does not rely on low-rank assumptions, the resulting GMM estimator is correspondingly more robust than those obtained using existing approaches.

It is technically challenging to characterize how imputation error affects the asymptotic distribution of the GMM estimator. This difficulty arises because (i) imputation error is non-classical and correlated across links, and (ii) the moment conditions rely on network statistics, which aggregate these errors and thereby further complicate the analysis. In this paper, we focus on the widely used linear-in-means peer-effects model and derive the asymptotic orders of the bias and variance of the resulting GMM estimator. We show that the optimal bandwidth for link-level imputation is generally not optimal for downstream estimation. We also provide guidance to help applied researchers assess the valid inference of parameter estimates obtained using our imputation approach.

We provide extensive simulation evidence to evaluate the finite-sample performance of our imputation method. First, we compare its imputation accuracy with alternative approaches and show that our method outperforms them across a wide range of sampling rates and network sparsity levels. We then evaluate the performance of our method in empirical analysis, including regressions based on network statistics and the linear-in-means peer-effects model. These simulations demonstrate that the estimates of regression coefficients and GMM estimators based on our imputed networks exhibit strong finite-sample performance. Finally, we illustrate the practical usefulness of the proposed method through an empirical application to the study of how social networks affect microfinance adoption banerjee2013diffusion.

\paragraph{Related literature} We contribute to the literature on estimating model parameters with sampled network data by developing an imputation method that allows for a fully nonparametric network formation model. chandrasekhar2011econometrics show that models estimated using sampled networks are generally biased. They propose an analytical correction for certain special cases and, for more general settings, they estimate a parametric network formation model and conduct GMM using simulated moments. boucher2020estimating study the estimation of the linear-in-means peer-effects model when the network data are corrupted, taking egocentrically sampled networks as a special case. Their approach relies on obtaining a consistent estimator of the distribution of the network, which is typically available only under parametric specifications. herstad2023estimating also examine peer-effects estimation with missing links, focusing on settings with large networks and relying on the network formation model of graham2017econometric. A key distinction relative to this line of work is that our method does not rely on parametric specifications of the network formation process and does not require specifying the distribution of unobserved heterogeneity. Therefore, parameter estimation based on our imputation method is more robust. hsieh2024non propose weighted estimators to recover several network-level statistics without assuming any network formation model. In comparison, our approach accommodates more general sampling mechanisms and supports a broader class of downstream regression analyses. In addition, our approach to GMM estimation differs from that of chandrasekhar2011econometrics, who construct moments by integrating over the missing network data. Instead, we directly plug the imputed network into the moment conditions to avoid heavy computational burdens.

Another related line of research studies the recovery of missing links under partial network data or even in the absence of network data. mccormick2015latent and breza2020using demonstrate that aggregated relational data can help recover network structure without observing any links. de2025identifying study how network ties can be identified from panel data. thirkettle2019identification show that network statistics are set identified under partially sampled networks, although their analysis relies on a strategic network formation model, whereas ours is based on a dyadic formation framework. Other sampling designs have also been studied—for example, induced subgraphs in chatterjee2015matrix and censored networks in griffith2022name.

Our work also relates to the growing literature on estimation with mismeasured network data, as empirical analysis based on imputed networks inevitably inherits the errors introduced in the imputation stage. For peer-effects models, lewbel2024ignoring show that when the number of misclassified links grows at a rate slower than the network size, the 2SLS estimator of bramoulle2009identification remains consistent and standard inference methods remain valid. lewbel2025estimating propose an adjusted 2SLS estimator for peer-effects models under i.i.d. link misclassification. hardy2019estimating also study link misclassification under i.i.d. errors and obtain consistent treatment effects using an expectation-maximization algorithm. cai2022linear examine linear regression on centrality measures under link-level i.i.d. Bernoulli errors, and propose bias-correction and inference methods for OLS estimators. The setting we consider is technically more challenging than those in the existing literature. Our imputed networks contain errors that are dependent across links and may have bias of the same order as their standard deviation. Our paper contributes to the literature by characterizing the convergence rate of linear-in-means peer-effects estimators in settings that allow for far more general forms of network measurement error.

Methodologically, the studies most directly connected to ours are in the recent literature on nonlinear factor models in panel-data causal inference, particularly the methods of feng2020causal, feng2023optimal, and deaner2025inferring, which also allow for nonparametric and nonlinear factor structures. Our local two-way fixed-effects regression is analogous to a local linear regression and therefore achieves a second-order approximation error. In contrast, the estimator in deaner2025inferring corresponds to a local constant regression and attains only a first-order approximation error. The local PCA methods of feng2020causal and feng2023optimal can be viewed as local polynomial regressions that fully exploit smoothness of nonlinear functions. Our approach differs from local PCA in that it does not require specifying the latent factor dimension. Another important difference is that, unlike most methods in this strand of the literature, our approach flexibly incorporates observed covariates to improve prediction accuracy.

Finally, our analysis is connected to several strands of the econometrics and statistics literature. The use of pseudo-distance as a measure of latent similarity was first developed by zhang2017estimating and has since been applied in a range of contexts, including nonparametric graphon estimation zhang2017estimating, zeleneev2020identification, controlling for unobservables using network data auerbach2022identification, and nonlinear factor models feng2020causal, feng2023optimal, beyhum2024inference, mugnier2025simple, deaner2025inferring. Our analysis is also related to the dyadic network formation literature graham2017econometric, dzemski2019empirical, chen2021nonlinear, ma2022detecting and the graphon estimation literature gao2015rate, xu2018rates. kitamura2024estimating consider covariate-assisted graphon estimation similar to ours, but they focus on stochastic block models, whereas we focus on more general nonlinear link-formation process. Finally, extracting the component of the network explained by covariates relies on the dyadic nonparametric regression developed by graham2021minimax.

The rest of the paper is organized as follows. Section (ref) introduces the framework and estimation procedure. Section (ref) provides statistical guarantees for the imputation. Section (ref) presents theoretical properties of the GMM estimators based on imputed networks. Section (ref) presents simulation results for both the imputation and the downstream empirical analysis. The empirical application to banerjee2013diffusion is given in Section (ref). Section (ref) concludes.

Method

Network Formation Model

Let $A\in\{0,1\}^{N\times N}$ denote the adjacency matrix of an unweighted network, where $A_{ij} = 1$ if there is a link between $i$ and $j$ and $A_{ij} = 0$ otherwise. We impose $A_{ij} = A_{ji}$ for all $i, j$ and $A_{ii} = 0$, so that the network is symmetric and has no self-loops. For each individual $i$, let $X_i \in \mathbb{R}^{d_X}$ denote the observed covariates (e.g., gender, education, and other demographic characteristics), and let $\xi_i \in \mathbb{R}^{d_\xi}$ denote the latent factors. The network $A$ is generated according to the following process. Conditional on observable and unobservable characteristics $\{X_i,\xi_i\}_{i=1}^{N}$, links $\{A_{ij}\}_{1\leq i<j\leq N}$ are generated independently with probability $\mathbb{P}\left(A_{ij} = 1 \mid \{X_i,\xi_i\}_{i=1}^{N} \right) = f(X_i, \xi_i, X_j, \xi_j)$, where $f: \mathbb{R}^{d_X + d_{\xi}}\times \mathbb{R}^{d_X + d_{\xi}} \rightarrow [0, 1] $ is a symmetric function, referred to as the graphon, whose functional form is unknown to econometricians. It is therefore convenient to write

where $\epsilon_{ij}$ are realization errors, which is conditionally independent across pairs $\{(i,j) \mid 1\leq i< j \leq N\}$ given $\{(X_i,\xi_i)\}_{i=1}^{N}$.

The network formation model in (ref) is a conditionally independent dyadic model. Conditional on the observable and unobservable characteristics $\{X_i,\xi_i\}_{i=1}^{N}$, the graphon $f$ does not depend on other individuals or on other links in the network. In other words, while $A_{ij}$ and $A_{ik}$ ($k \neq j$) may be correlated unconditionally, they are independent given $\{X_i,\xi_i\}_{i=1}^{N}$. This framework is flexible and encompasses a wide range of commonly used network formation models. Below, we provide several examples.

Beyond these examples, our network formation framework in (ref) can flexibly accommodate richer forms of heterogeneity. For example, it allows for heterogeneous slope coefficients when the homophily effect on observables, captured by $\beta$ in (ref), depends on the latent factors $(\xi_i, \xi_j)$. Similarly, the unobserved heterogeneity term $g(\xi_i, \xi_j)$ can capture more complex interaction structures, such as interactive fixed effects of the form $\xi_i'\xi_j$. However, our dyadic model excludes strategic network formation, where the formation of one link may depend on the presence of other links (see de2020econometric).

Data

Economists observe individual covariates $X_i$ (e.g., gender, age) for all $i$ in the network but do not observe the full network $A$; instead, only a sampled network is available. This setting is common in applications because collecting covariates is typically less costly than collecting network information.

In this paper, we focus on egocentric sampling, which is one of the most commonly used sampling designs in empirical work. Specifically, we randomly sample $n < N$ nodes in the network and ask all their friends not only in the sample but in the whole network. The sampling probabilities can depend on both covariates $X_i$ and unobserved characteristics $\xi_i$. Let $\mathcal{S} \subset \{1,2, \ldots, N\}$ denote the set of sampled individuals with $|\mathcal{S}| = n$. Then a link $A_{ij}$ is observed if $i\in \mathcal{S}$ or $j\in \mathcal{S}$, whereas links between two non-sampled nodes (i.e., $i\notin\mathcal{S}$ and $j\notin\mathcal{S}$) are unobserved. A simple example follows to illustrate the setup.

Since it causes no harm to relabel the agents, we can do so such that $\mathcal{S} = \{1,2,\ldots,n\}$, in which case the observed adjacency matrix $A^{\mathrm{obs}}$ admits the following block structure:

Here, $A_{\mathcal{S}\mathcal{S}}$ consists of links between two sampled individuals, $A_{\mathcal{S}^c\mathcal{S}} = A_{\mathcal{S}\mathcal{S}^c}' $ contains links between sampled and unsampled individuals, and the bottom right block $A_{\mathcal{S}^c\mathcal{S}^c}$ is unobserved because both individuals in each dyad are unsampled; we highlight this block in red.

Imputation

To illustrate the main idea of our imputation approach, let us start with a simple case in which both $X_i$ and $\xi_i$ are observable for every individual in the network. In this case, imputing missing links would be straightforward: we could estimate $f(\cdot, \cdot)$ nonparametrically from the observed adjacency matrix to obtain estimates $\hat{f}(\cdot, \cdot)$, and then use these estimates to impute the missing links in the submatrix $A_{\mathcal{S}^c \mathcal{S}^c}$.

The main difficulty, however, is that the latent factors $\xi$ are unobserved for all individuals. To make this challenge more transparent, let us first abstract from the observed covariates $X$ and focus on the case where link formation depends solely on the latent factors. In this setting, the imputation problem reduces to recovering information about the latent factors $\xi$ from the observed subnetwork and using them to predict the missing entries of the adjacency matrix $A_{\mathcal{S}^c \mathcal{S}^c}$.

Our approach is motivated by the following idea. When network formation is dyadic, as in (ref), individuals with similar $\xi$ are expected to have similar linking probabilities. This observation provides a useful heuristic: similarity in observed connections reflect similarity in latent positions. Building on this intuition, for any unsampled individual $i$, we use their observed links to sampled individuals to measure their similarity to sampled nodes. We then use the observed links of sampled individuals to other unsampled individuals, weighted by their similarity to individual $i$, to impute $i$'s missing links to other unsampled individuals.

Because the information about $\xi$ is estimated rather than directly observed, it is inappropriate to ignore the observed covariates or absorb them into the latent factors, as doing so would discard the noiseless information contained in $X$. Therefore, when $X$ is available, our method first extracts the component of link formation explained by the observed covariates and then applies the latent-factor-based procedure described above to the residuals. These two components are subsequently combined to obtain the final imputed network.

Extract information explained by $X$

As discussed previously, when individual covariates $X$ are available, it is appropriate to first extract the component that can be explained by $X$. Consider the decomposition:

where $\Pi_{ij}$ represents the part of the graphon that can be explained by the observed covariates, and $\Pi_{ij}^{\bot}$ captures the residual component of the graphon that cannot be explained solely by the observed covariates. Note that the observed link $A_{ij}$ can then be written as $A_{ij} = \Pi_{ij} + (\Pi_{ij}^{\bot} + \epsilon_{ij})$, and the composite error term satisfies $\mathbb{E}(\Pi_{ij}^{\bot} + \epsilon_{ij}\mid X_i, X_j) = 0$, which implies that dyadic nonparametric regression of $A_{ij}$ on $(X_i, X_j)$ consistently estimates $\Pi_{ij}$. graham2021minimax derive the pointwise and uniform minimax risks for estimating such conditional expectations and show that the Nadaraya-Watson (NW) estimator attains the optimal convergence rate\footnote{ In practice, we implement local linear regression rather than using NW estimator. We discuss this in Appendix (ref). }. We implement the dyadic nonparametric regression using the observed subset of the adjacency matrix. Since the covariates $X$ are also observed for unsampled individuals, the nonparametric estimator can then be evaluated directly for missing links, that is, the estimator $\hat{\Pi}_{ij}$ can be computed even when both $i$ and $j$ are unsampled. The implementation details are discussed in Appendix (ref).

In the next two steps, we focus on estimating this residual component ${\Pi}^{\bot}$ using pseudo-distance and two-way fixed-effect imputation procedure described below.

Calculate similarity

We now describe how to infer similarity between latent factors using observed network links.

For notational simplicity, let $\zeta_{i} := (X_{i}, \xi_{i})$ for each $i$, and let $f_{ij}=f(\zeta_i, \zeta_j)$ for all $1\leq i, j\leq N$. The adjacency matrix $A$ can be expressed as $A = P + E$. Here, $P$ collects the conditional probabilities of link formation, i.e., $P_{ij} = f_{ij}$ when $i\neq j$, and $P_{ii} = 0$ (because self-loops are ruled out), and $E$ collects error terms $\epsilon_{ij}$\footnote{ Likewise, we impose $\epsilon_{ii} = 0$ for each $i$. }. Our approach relies on the following idea. For any unsampled individual $i$ and sampled individual $i'$, if their linking probabilities to other sampled individuals (or equivalently, the corresponding rows in $P$) appear sufficiently similar, then under suitable regularity conditions, $\zeta_i$ and $\zeta_{i'}$ are also similar in the latent space.

A natural measure of the difference between linking probabilities is the squared $L_2$ distance between the functions $f(\zeta_i, \cdot)$ and $f(\zeta_{i'}, \cdot)$, defined as

The squared $L_2$ distance equals zero if and only if the two individuals $i$ and $i'$ have identical link formation probabilities with all other individuals in the network, that is, $f(\zeta_i, \cdot) = f(\zeta_{i'}, \cdot)$ almost everywhere. Therefore, if one could construct a consistent estimator of $L^2_2(\zeta_i, \zeta_{i'}) $ from the observed network, it would provide a natural proxy for latent similarity. However, the squared $L_2$ distance is not directly estimable. Its sample analogue based on observed network links is given by

where the additional variance terms do not vanish asymptotically and therefore contaminate the estimation. A formal derivation and further discussion are provided in Appendix (ref).

To address this issue, we follow zhang2017estimating and define a new metric between the functions $f(\zeta_i, \cdot)$ and $f(\zeta_{i'}, \cdot)$ as

The new metric $d(\cdot, \cdot)$ has two useful properties. First, it is straightforward to verify that $d(\zeta_i, \zeta_{i'}) \geq \frac{1}{2} L^2_2(\zeta_i, \zeta_{i'})$ and that $d(\zeta_i, \zeta_{i'}) = 0$ if $\zeta_i = \zeta_{i'}$. Therefore, controlling $d(\cdot,\cdot)$ is sufficient to control $L^2_2(\zeta_i, \zeta_{i'})$. Second, unlike the squared $L_2$ distance, $d(\zeta_i, \zeta_{i'})$ admits a consistent estimator based on the observed network. Specifically, define

We show that $\hat{d}_{ii'} \stackrel{p}{\longrightarrow} d(\zeta_i, \zeta_{i'})$. We refer to $\hat{d}_{ii'}$ as the pseudo-distance and $d(\zeta_i, \zeta_{i}')$ as the population pseudo-distance throughout the paper. Calculating $\hat{d}_{ii'}$ relies only on the friendships of sampled individuals and can therefore be computed from the observed subnetwork. In addition, under suitable regularity conditions, the pseudo-distance $\hat{d}_{ii'}$ is informative about latent similarity, in the sense that small values of $\hat{d}_{ii'}$ imply that $\zeta_i$ and $\zeta_{i'}$ are close.

It is important to note that the discussion here concerns the composite characteristic $\zeta_i$, rather than the unobserved heterogeneity $\xi_i$ alone. This is because (i) $X_i$ and $\xi_i$ enter the network formation process in a non-separable way, making it impossible to isolate their effects without additional assumptions; and (ii) the information in $X_i$ that can be directly explained has already been extracted in the dyadic nonparametric regression, so there is no need to distinguish between $X_i$ and $\xi_i$ in this step.

Local two-way fixed-effects imputation

For the projection residuals, we employ the pseudo-distance $\hat{d}$ obtained in the previous step to construct a local estimator for $\Pi_{ij}^{\bot}$, referred to as the two-way fixed-effect imputation. Specifically, for unsampled individuals $i$ and $j$, we estimate the local two-way fixed-effects regression to obtain $\hat{a}_i$ and $\hat{b}_j$:

where $K_h(\cdot) := h^{-1}K(\cdot / h)$ is a kernel function with bandwidth $h$, and $\hat{d}$ denotes the pseudo-distance. Then the imputed $\Pi_{ij}^{\bot}$ is obtained via $\hat{\Pi}_{ij}^{\bot} = \hat{a}_i + \hat{b}_j$, and the imputed network link is given by

The imputation is “local” in the sense that, by using a kernel estimator, we rely only on information from individuals $i'$ and $j'$ whose pseudo-distances $\hat{d}_{ii'}$ and $\hat{d}_{jj'}$ (or, equivalently, the distances between their latent positions) are sufficiently small when estimating fixed effects. Note that, although our notations omit the indices $(i, j)$ for simplicity, the estimator $(\hat{\boldsymbol{a}}, \hat{\boldsymbol{b}})$ in (ref) is pair-specific and therefore must be re-estimated for each $(i, j)$.

The procedure is referred to as a “two-way fixed-effects imputation" because (ref) takes the form of the two-way fixed-effects regression rather than a local constant regression. This is analogous to replacing a local constant estimator with a local polynomial estimator in nonparametric regression: incorporating the two-way fixed-effects provides a higher-order local approximation and thus improves estimation accuracy. We will formally justify this result later using a Taylor expansion.

It is important to note that $\hat{A}_{ij}$ can only approximate the conditional link-formation probability $f_{ij}$ rather than the missing link $A_{ij}$. This is because the missing link contains an unpredictable idiosyncratic error $\epsilon_{ij}$, which cannot be predicted from the available data.

We impute all missing links based on the previous steps and merge these imputed entries with the observed ones to obtain a complete imputed matrix $\hat{A}$. To facilitate understanding of the main idea of our algorithm, we summarize below a simplified version of our method that excludes cross-fitting and cross-validation, which we would discuss later.

Cross-validation and sample-splitting

Although the proposed algorithm is straightforward to implement, deriving the asymptotic properties of $\hat{A}_{ij}$ can be technically challenging because of the potential dependence between the estimated pseudo-distance $\hat{d}$ and the error terms $\epsilon$ in the two-way fixed-effects regression.

A standard approach to eliminate this dependence is to employ cross-fitting, in which the data are partitioned into folds so that the first-step and second-step estimations are conducted on separate subsets of the data. Specifically, we randomly split the sampled individuals $\mathcal{S}$ into two parts, $\mathcal{S}_1$ and $\mathcal{S}_2$, of equal size. We then use only the friendships of $\mathcal{S}_1$ to calculate the similarity $\hat{d}(\cdot,\cdot)$:

In the imputation stage, we use no link information from $\mathcal{S}_1$, and all calculations rely on the friendships of $\mathcal{S}_2$:

In addition, a practical concern is the choice of the tuning parameter $h$. We recommend selecting $h$ using a leave-one-out cross-validation procedure, and the details of cross-validation, together with the sample-splitting, are provided in Appendix (ref).

Theory for Imputation

We now introduce the regularity conditions.

Assumption (ref)(ref) indicates that our method is designed for large-network settings in which both the network size $N$ and the number of sampled individuals $n$ grow large. Our analysis does not rely on $n$ and $N$ growing at the same rate and continues to hold when $n$ grows slowly, e.g., $n = O(N^{\alpha})$ for any $\alpha \in (0,1)$.

Assumption (ref)(ref) imposes that $\{(X_i, \xi_i)\}_{i=1}^{N}$ are i.i.d., and the data-generating process of $A$ follows a conditionally independent dyadic model (ref). This model also follows from the Aldous-Hoover Theorem aldous1981representations, hoover1979relations, which shows that any infinitely exchangeable random graph admits such a representation. As discussed earlier, although (ref) is flexible enough to encompass a wide range of network formation models commonly used in economics and statistics, it rules out strategic network formation, in which individuals linking decisions depend on those of others; see de2020econometric for a review.

Assumption (ref)(ref) is standard in egocentric sampling designs, which assume that individuals are randomly drawn from the population. Our method also applies to more general sampling schemes in which individuals are sampled with heterogeneous probabilities that depend on both observed and unobserved characteristics, i.e., $D_i \bot A \mid \{(X_i, \xi_i)\}_{i=1}^{N}$. In this case, the asymptotic analysis can be extended with only minor modifications, as long as the sampling probabilities $\mathbb{P}(D_i = 1\mid \{(X_i, \xi_i)\}_{i=1}^{N})$ are uniformly bounded away from zero across individuals.

The first part of Assumption (ref)(ref) imposes the boundedness of $\zeta$, and the second part introduces a standard bounded-density condition. Intuitively, it ensures that the distribution of $\zeta_i$ is neither locally too sparse nor too concentrated by requiring that, for any point $\zeta$ and a small radius $h>0$, the probability that $\zeta_i$ falls within a ball of radius $h$ centered at $\zeta$ is of the same order as the volume of that ball. This requirement is automatically satisfied when (i) $\mathrm{supp}(\zeta)$ is compact and convex, (ii) $\mathrm{affine}(\mathrm{supp}(\zeta)) = \mathbb{R}^{d_{\zeta}}$, and (iii) the probability density function of $\zeta$ is uniformly bounded above and bounded away from zero. This assumption is common for establishing the asymptotic properties of nonparametric kernel estimators.

Assumption (ref)(ref) requires the graphon $f$ to be twice continuously differentiable with uniformly bounded second derivatives. This smoothness condition is needed because our two-way fixed-effects imputation method relies on a second-order Taylor expansion of $\Pi^{\bot}_{ij}$ within a local neighborhood to eliminate the first-order approximation bias, thereby improving accuracy. Methods that rely on higher-order Taylor expansions, such as the local PCA approach proposed by feng2023optimal, correspondingly require the existence of higher-order derivatives. We discuss this point in detail later.

It should be noted that the requirement of second-order differentiability of $f$ is not a mild condition. For instance, consider the network formation model with degree heterogeneity in Example (ref). If the distance function $\omega(\cdot,\cdot)$ is defined by the absolute distance, i.e., $\omega(X_i, X_j) = |X_i - X_j|$, which is common in applied research, then the corresponding graphon $f$ is at most Lipschitz continuous rather than twice differentiable. Nevertheless, our method remains valid and can consistently impute the conditional probability of link formation, although the approximation error no longer attains the second-order rate.

Asymptotic analysis of pseudo-distance

We now develop the asymptotic properties of the pseudo-distance $\hat{d}_{ii'}$ and show that a sufficiently small $\hat{d}_{ii'}$ implies that $\zeta_i$ and $\zeta_{i'}$ are close in the latent space. In the following lemma, we establish the convergence of $\hat{d}_{ii'}$ to its population counterpart $d(\zeta_i, \zeta_{i'})$.

Lemma (ref) shows that the approximation error of $\hat{d}_{ii'}$ relative to its population counterpart $d(\zeta_i, \zeta_{i'})$ can be decomposed into two components. The first term, $\gamma_1 \frac{\log N}{N^{1/d_{\zeta}}}$, captures the matching discrepancy in the latent characteristics $\zeta_i$. Its rate deteriorates with the dimension of the latent space $d_{\zeta}$, which is consistent with the well-known “curse of dimensionality” in the matching literature. The second term, $\gamma_2 \frac{\log N}{\sqrt{n}}$, reflects the effect of the idiosyncratic error $\epsilon_{ij}$. It shows that although the true graphon is unobservable, the pseudo-distance can asymptotically denoise the observed network. The matching discrepancy error dominates the effect of idiosyncratic errors when $N\sim n$ and $d_{\xi} \geq 3$, which is the common case in practice. Finally, for notational simplicity, let $$\delta_{N, n} := \gamma_1 \frac{\log N }{N^{1/d_{\zeta}}} + \gamma_2 \frac{\log N }{\sqrt{n}}, $$ and focus on the event under which $|\hat{d}_{ii'} - d(\zeta_i, \zeta_{i'})|$ is uniformly bounded by $\delta_{N, n}$, as implied by Lemma (ref).

The following assumption is critical and requires that the population pseudo-distance $d(\zeta_i, \zeta_{i'})$ be sufficiently informative to reveal the latent distance $\|\zeta_i - \zeta_{i'}\|$.

This assumption ensures that when the population pseudo-distance $d(\zeta_i, \zeta_{i'})$ is sufficiently small, the latent positions $\zeta_i$ and $\zeta_{i'}$ are close to each other. Hence, $d(\zeta_i, \zeta_{i'})$ can be used to identify individuals with similar latent characteristics. This assumption is crucial for establishing link-level imputation error bounds. In Appendix (ref), we provide a set of primitive sufficient conditions for Assumption (ref) and verify that it holds under commonly used network formation models.

\setcounter{example}{0}

\paragraph{Remark} An alternative to the pseudo-distance is the difference in average degrees, defined as $\hat{\rho}_{ii'} = \frac{1}{|\mathcal{S}_1|}\sum_{\ell \in \mathcal{S}_1}A_{i\ell} - \frac{1}{|\mathcal{S}_1|}\sum_{\ell \in \mathcal{S}_1} A_{i'\ell}$. However, using the average-degree distance requires a stronger informativeness assumption than Assumption (ref). In particular, even when $\zeta$ is one-dimensional, using average-degree distance requires that the population average degree, $\zeta \mapsto \int f(\zeta, \tilde{\zeta})\,d\mathbb{P}(\tilde{\zeta})$, is injective in $\zeta$. (In the case where $\zeta$ is a scalar, this is equivalent to requiring the population average degree to be strictly monotonic in $\zeta$.) This condition does not hold in either Example (ref) or Example (ref) and is stronger than Assumption (ref) in our paper.

\paragraph{Remark} While Assumption (ref) is convenient for theoretical analysis, it should not be interpreted as necessary for good imputation performance. Even under the cases where the informativeness condition may fail, our simulation results in Appendix (ref) show that our proposed method still performs well in practice and outperforms alternative approaches.

Asymptotic analysis of imputation

The error of our proposed imputation method can be decomposed into two parts corresponding to the dyadic nonparametric regression step and the local two-way fixed-effects regression step. The first-stage error is typically asymptotically negligible relative to the second-stage error for two reasons. First, the dyadic nonparametric regression only extracts variation explained by the observed covariates $X$, whereas the local two-way fixed-effects regression involves unobserved heterogeneity $\xi$. Because $\xi$ must be approximated through the pseudo-distance, there is an additional source of estimation error in the second stage. Second, the dimension of the regressors in the first-stage nonparametric regression is $2d_X$. By contrast, because the graphon $f$ is nonseparable in $X$ and $\xi$, the dimension in the second-stage regression is $d_{\zeta} = d_X + d_{\xi}$, leading to a slower convergence rate. Therefore, our asymptotic analysis primarily focuses on the second-stage estimation error.

The asymptotic properties of dyadic nonparametric regression are studied in the literature (see graham2021minimax). We therefore summarize the required conditions in the following assumption.

Assumption (ref) imposes uniform convergence rates on the bias and variance of the first-stage dyadic nonparametric regression estimator. The convergence rates depend only on $d_X$ and the number of sampled individuals $n$ (up to a logarithmic factor). Under Assumption (ref), these rates can be achieved by the Nadaraya-Watson estimator, as shown in graham2021minimax.

\paragraph{Remark} The convergence rates in Assumption (ref) are the MSE-optimal rates for twice continuously differentiable target functions with $n$ i.i.d. observations. As noted by graham2021minimax, dyadic nonparametric regression differs from standard nonparametric regression with i.i.d. data in two respects: (i) although the dyadic conditional expectation $\Pi(\cdot,\cdot)$ takes $2d_X$ arguments, the convergence rates behave as if $\Pi$ depended on only $d_X$ arguments; and (ii) when the entire network is observed, although the regression involves $N^2$ dyads, the effective sample size is just $N$ because the estimator is a U-statistic. The convergence rates in Assumption (ref) are consistent with the first feature, but differ from the second because the network is incomplete in our setting. One use Hoeffding decomposition to show that the effective sample size in our setting is the number of sampled individuals $n$ rather than $N$.

This assumption imposes regularity conditions on kernel $K$ used in the second step. The requirement is standard in the literature, and many commonly used kernels, including the Epanechnikov kernel, satisfy Assumption (ref).

We now briefly explain why $\hat{d}_{ii'}$ can be viewed as a noisy measure of $\|\zeta_i - \zeta_{i'}\|$ and can thus serve as its proxy in the kernel function. First, let $h$ denote the bandwidth of the kernel function. When imputing $\Pi^{\bot}_{ij}$, Assumption (ref) allows us to restrict attention to links $(i', j')$ for which both $\hat{d}_{ii'}\leq h$ and $\hat{d}_{jj'}\leq h$. Second, by Lemma (ref), we have $|\hat{d}_{ii'} - d(\zeta_i, \zeta_{i'})| \leq \delta_{N, n}$ with high probability. When $\delta_{N, n} / h \to 0$, the approximation error in $\hat{d}_{ii'}$ is asymptotically negligible relative to the bandwidth $h$, so $d(\zeta_i, \zeta_{i'})$ can be replaced by $\hat{d}_{ii'}$ without loss of first-order accuracy. Furthermore, under Assumption (ref), $d(\zeta_i, \zeta_{i'})$ is informative about the latent distance in the sense that $\|\zeta_i - \zeta_{i'}\|\leq \lambda^{-1} d(\zeta_i, \zeta_{i'})$. Combining these results yields, with probability approaching one, $\|\zeta_i - \zeta_{i'}\|\leq \lambda^{-1} h$, which implies that using pseudo-distances in the kernel function selects pairs whose latent distances lie within a $\lambda^{-1} h$-neighborhood when imputing missing links.

The following theorem provides the theoretical guarantee for our imputation method.

Theorem (ref) establishes link-level error bounds, rather than a Frobenius-norm bound for a submatrix (e.g., gao2015rate and xu2018rates). The convergence is uniform over all missing links, which implies that the imputation error is controlled at the link level and is particularly useful for applications that require accurate imputation of each link probability.

Following the classical nonparametric regression literature, we decompose the estimation error into bias and variance. The bias term is of order $\left(h^2 + 1/(n h^{d_{\zeta}})\right)$, where (i) $h^2$ arises because our local two-way fixed-effects regression removes the first-order bias, leaving only the second-order approximation error, and (ii) $1/(n h^{d_{\zeta}})$ reflects the contamination from diagonal entries of the matrix $A$ (which are $0$ instead of $f_{ii} + \epsilon_{ii}$ because we rule out self-loops). In most empirically relevant cases, the first term, $h^2$, is dominant. In particular, under the MSE-optimal bandwidth choice, it can be directly verified that $1/(n h^{d_{\zeta}}) = o(h^2) $.

The variance consists of two terms. The first term, $\log N /(n h^{d_{\zeta}})$, matches the rate in the standard nonparametric literature. The second term, $\delta^2_{N,n} h^2$, is specific to our estimator and arises from the use of the pseudo-distance as a substitute for the unobserved latent factors $\zeta_i$ and $\zeta_j$. When the latent factors are observable, we have $\delta_{N,n}=0$, and the second term disappears, reducing our result to classic conclusions in the nonparametric regression literature (see Tsybakov2008introduction and, for dyadic regression, graham2021minimax).

The bandwidth $h$ governs the bias-variance trade-off. However, because the bias and variance each consist of two components, the trade-off deviates from the classical setting. As $h$ decreases, $h^2$ and $h^2 \delta_{N, n}^2$ diminish, while $1/(n h^{d_{\zeta}})$ increases. Therefore, the optimal choice of $h$ depends jointly on the relative magnitudes of $n$ and $N$ and on $d_{\zeta}$. Nevertheless, it is straightforward to verify that, when $n$ and $N$ are of the same order, setting $h \asymp n^{-1/(4 + d_{\zeta})}$ attains the optimal MSE rate (up to logarithmic factors).

\paragraph{Remark} Our method employs two-way fixed-effects regression for imputation, which effectively removes the first-order approximation error. The intuition behind this approximation can be illustrated by a Taylor expansion. Heuristically, to impute $A_{ij}$, if we ignore the first-step nonparametric dyadic regression and consider the Taylor expansion of $f(\zeta_{i'}, \zeta_{j'})$ around the fixed point $(\zeta_i, \zeta_j)$:

Because we employ a kernel-weighted estimator based on the pseudo-distances $\hat{d}_{ii'}, \hat{d}_{jj'}$ with bandwidth $h$, and since, as discussed above, $\hat{d}_{ii'}$ and $\hat{d}_{jj'}$ differ from the latent-space distances $\|\zeta_{i'} - \zeta_{i}\|$ and $\|\zeta_{j'} - \zeta_{j}\|$ only up to a bounded factor, the contribution to imputation comes from links $(i', j')$ with $\hat{d}_{ii'}, \hat{d}_{jj'}\leq h$. For those links, we obtain

Hence, the two-way fixed-effects regression absorbs first-order approximation errors, leaving only the second-order remainder $O(h^2)$.

The idea of using two-way fixed-effects imputation is closely related to the approaches developed in freeman2023linear and beyhum2024inference. Our approach differs from theirs in two respects. First, their methods do not involve missing data and are primarily designed to estimate structural parameters in panel data models, whereas our method aims to predict the formation probability of missing links. Second, they discretize the data into groups using k-means clustering, while we employ a kernel estimator.

The methods most conceptually related to ours are those of feng2023optimal and deaner2025inferring. deaner2025inferring propose a local-constant estimator to impute potential outcomes. Because they assume only that the potential outcome function is Lipschitz continuous, their approximation error is of order $O(h)$. Building on the ideas of freeman2023linear and beyhum2024inference, we assume that the graphon is twice continuously differentiable and, by using a two-way fixed-effects imputation, we reduce the approximation error to $O(h^2)$. In comparison, the local PCA method proposed by feng2023optimal leverages factor models together with higher-order Taylor expansions to achieve more accurate approximations. In particular, feng2023optimal show that the local PCA estimator attains Stone's minimax-optimal rate for nonparametric regression. We refer interested readers to that paper for further details.

We adopt the two-way fixed-effects model instead of local PCA approach for two reasons. First, local PCA method requires prior knowledge of the dimension $d_{\zeta}$, which is typically unknown in practice. In contrast, our approach does not require prior knowledge of $d_{\zeta}$. Second, networks in most economic applications are relatively small. For instance, in development economics, village networks typically contain around $200$ households, and with a sampling rate of roughly $40\%$, the effective sample size available for imputing a specific missing link can be fewer than $20$. Under such small-sample conditions, achieving theoretical higher-order accuracy is difficult. Our numerical simulations also show that the finite-sample performance of our method is better to that of the local PCA.

Theory for Downstream Estimation

We consider the application of our method to estimating parameters in downstream economic models. The data $\{(A_m, X_m, W_m, Y_m)\}_{m=1}^{M}$ consist of $M$ large i.i.d. networks indexed by $m = 1,2,\ldots, M$, and each network has size $N_m$. Here, $A_m\in \{0, 1\}^{N_m\times N_m}$ is the adjacency matrix, $Y_m \in \mathbb{R}^{N_m}$ is the outcome variable, $X_m \in \mathbb{R}^{N_m\times d_{X}}$ denotes the covariates that directly enter the network formation model as in (ref), whereas $W_m \in \mathbb{R}^{N_m\times d_{W}}$ contains variables that do not necessarily appear in the network formation model but affect the outcome $Y_m$. Note that $W_m$ and $X_m$ may overlap, i.e., $W_m\bigcap X_m \neq\emptyset$.

The adjacency matrix $A_m$ is not fully observed because of egocentric sampling, but the outcome $Y_m$ and the covariates $(X_m, W_m)$ are observed for all individuals in the network. Let $n_m$ be the number of sampled individuals in network $m$. We mainly focus on the asymptotics in which $M, N_m, n_m\rightarrow \infty$. The large-network asymptotics $N_m, n_m \rightarrow \infty$ are required for consistent imputation of the adjacency matrix. The large-$M$ asymptotics are introduced to accommodate within-network dependence of error terms in the downstream regression. If one is willing to impose independence of these errors within a network (see Examples (ref) and (ref) in the following text), then the consistency analysis can also be extended to the single-large-network setting.

We allow for heterogeneous network formation process $f_m$ for each $m = 1, \ldots, M$, so that information from other networks does not help in imputing a given network. Therefore, we implement imputation method separately for each incomplete network, obtain the imputed adjacency matrix $\hat{A}_m$, and use it, rather than the observed network $A^{\mathrm{obs}}$, in the moment conditions to contruct the GMM estimator.

Consistency

\paragraph{GMM} We work with the moment function $$\psi\left(A_m, W_m, Y_m, \alpha \right): [0, 1]^{N_m\times N_m} \times \mathbb{R}^{N_m\times d_{W}} \times \mathbb{R}^{d_{\alpha}} \mapsto \mathbb{R}^{q}, $$ where $q$ is the number of moment conditions and $\alpha\in \mathbb{R}^{d_{\alpha}}$ is the structural parameter of interest in the downstream analysis. We assume the following moment condition

where $\alpha_0$ is the true parameter. Since $A_m$ is partially observed due to sampling, we impute each sampled network separately to obtain $\{\hat{A}_m\}_{ m=1, \ldots, M}$, and then use the imputed networks to solve

where $\Sigma_{M} \in \mathbb{R}^{q\times q}$ is the weight matrix.

\paragraph{Remark} Our approach for GMM estimation differs from that of chandrasekhar2011econometrics. Specifically, chandrasekhar2011econometrics construct moments by integrating over the missing network data, using simulation to approximate (conditional) expectations. As a consequence, their implementation is computationally demanding in GMM settings and relies on parametric assumptions about the distribution of the error terms (see the discussion in chandrasekhar2011econometrics and boucher2020estimating). By contrast, our approach directly plugs the imputed network into the moment function $\psi(\cdot)$, avoiding simulation and distributional assumptions on the errors. A key advantage of chandrasekhar2011econometrics's approach is that it can accommodate complex network statistics, such as average path length, that cannot be computed using a simple plug-in strategy. However, even in such cases, our imputed adjacency matrix can still serve as the basis for the simulation-based GMM procedure, and their analytical framework (developed in Section 4.2 of chandrasekhar2011econometrics) can then be used to establish consistency of the downstream estimator.

We now introduce a set of regularity conditions on the sampling process and the moment function $\psi(\cdot)$ used to establish consistency of the downstream estimator.

Assumptions (ref)(ref)--(ref) are standard conditions to ensure the consistency of M-estimator. Assumption (ref) considers the asymptotic regime in which the number of networks $M$ grows large. Assumption (ref) requires each network to be independent and identically distributed across $m$. Assumption (ref)(ref) is an identification condition and ensures that the moment condition does not hold at $\alpha\neq\alpha_0$. Assumption (ref) requires that the weighting matrix $\Sigma_M$ converges to a positive definite matrix as $M\rightarrow\infty$. Assumption (ref)(ref) requires the parameter space to be compact, the moment function to vary continuously with $\alpha$, and the second moments of $\psi$ to remain uniformly bounded over the parameter space. Assumption (ref)(ref) imposes a Lipschitz condition requiring that perturbations in $P_m$ lead to changes in $\psi$ up to a factor $L(W_m, Y_m)$ with a finite second moment, uniformly over $\alpha$. Since our parameter of interest is $\alpha$ and $P_m$ serves as a nuisance parameter, this assumption can be interpreted as a smoothness condition on the nuisance parameters within the framework of semi-parametric M-estimation.

Assumption (ref)(ref) requires that replacing the adjacency matrix $A$ with the probability matrix $P$ in the moment function becomes asymptotically negligible as the network size $N$ grows. Unlike the previous assumptions which are standard in the M-estimation and GMM literature, Assumption (ref)(ref) is specific to our setting and is the key condition for establishing consistency of the plug-in estimator. This is because our imputation procedure can only recover the conditional link-formation probabilities, rather than the realized missing links themselves. Intuitively, the assumption requires that the idiosyncratic link-formation shocks average out asymptotically in the moment conditions. This condition is satisfied in many empirical applications based on network statistics, such as linear-in-means peer-effects models and regressions on network centralities.

\setcounter{example}{2}

Theorem (ref) shows that, in the downstream regression, the plug-in estimator $\hat{\alpha}$ based on our imputed networks is consistent. It is particularly relevant in empirical applications, where interest typically focuses on the consistency of parameter estimates constructed from the imputed network rather than on the accuracy of the imputed links themselves. The key advantage of our approach is that the consistency of $\hat{\alpha}$ does not depend on (i) specification of the network formation model and (ii) the dimension of latent heterogeneity. In addition, even when a researcher prefers to impose a specific parametric network formation model for computational simplicity or because of strong prior knowledge, our method still provides a benchmark for assessing the extent to which downstream estimators depend on parametric assumptions about the network formation process.

\paragraph{Remark} The analysis above focuses on large-$N$, large-$M$ asymptotics. However, under certain conditions, the downstream regression analysis can also be extended to the single-large-network setting (i.e., $M=1$). For example, in Examples (ref) and (ref), consistency of the downstream GMM estimator continues to hold when the error terms $\{e_{i}\}_{i=1,\ldots,N}$ are uncorrelated across individuals.

Analysis of linear-in-means peer-effect model

Although Theorem (ref) establishes the consistency of the plug-in estimator based on our imputed adjacency matrices, analyzing its convergence rate and inference is challenging. First, as discussed in cai2022linear, even when each link contains only i.i.d. classical measurement error, network statistics used in moment conditions typically aggregate errors across links, thereby inducing non-classical measurement error that complicates the analysis. Second, the measurement error generated by our imputation procedure is correlated across links, and because the imputation is nonparametric, its bias and variance may be of the same order, further complicating the analysis of convergence.

Since network links typically enter the moment conditions through network statistics, the asymptotic properties of imputed links established in Theorem (ref) could be different from those of the GMM estimators that rely on the imputed network. Therefore, the bandwidth optimal for imputing links may not be optimal for downstream analysis. In the following text, we focus on the linear-in-means peer-effects regression and study the asymptotic behavior of the corresponding GMM estimator. We also provide guidance for downstream regression analysis using our imputed networks.

Consider the linear-in-means peer-effects model (ref) in Example (ref). The weight matrix is $\Sigma_M$. Let $\hat{V}_m := [\boldsymbol{1}, \hat{G}_m Y_m, W_m, \hat{G}_m W_m]$, $\hat{Z}_m := [\boldsymbol{1}, W_m, \hat{G}_m W_m, \hat{G}^2_m W_m]$, and $\alpha := (\alpha_C, \alpha_{\bar{Y}}, \alpha'_{W}, \alpha'_{\bar{W}})'$. The GMM estimator is given by

Assumption (ref)(ref) requires the networks to be independent and identically distributed across $m$. However, we allow the error terms $e_{mi}$ to be correlated across individuals within the same network. Assumption (ref)(ref) imposes the strict exogeneity condition and rules out correlated effects (bramoulle2020peer). Note that Assumption (ref)(ref) is stronger than the standard assumption in the literature, i.e., $\mathbb{E}\left(e_m \mid G_m, W_m\right) = 0$. This is because the network data are incomplete in our setting, requiring the error term in the downstream regression to be orthogonal to the information used for network imputation. We further assume that the error terms in the linear-in-means peer-effects model, $e_m$, are conditionally independent of the errors in the network formation process, $\epsilon_m$, which ensures that the imputation error is independent of $Y_m$.

Assumption (ref)(ref) requires that all networks have sizes of the same order and that the sample sizes within networks are also of the same order. In other words, no network or sample becomes disproportionately large or disproportionately small.

Assumption (ref)(ref) imposes a non-degeneracy condition requiring that the smallest singular value of $\mathbb{E}\left(\frac{1}{N_m} Z_m'V_m \right)' \mathbb{E}\left(\frac{1}{N_m} Z_m'V_m \right)$ is bounded away from zero. This requirement is analogous to the rank condition in instrumental-variable estimation, and is particularly important under large-network asymptotics. For example, if $W$ is i.i.d. within each network and independent of the network (e.g., random treatment assignment), the instrument $Z_m$ becomes asymptotically collinear, causing the smallest singular value of the above matrix to shrink toward zero and leading to an inconsistent IV estimator. In practice, one can let $W_m$ to be correlated with $X_m$ to satisfy this assumption (see hayes2024peer for further discussion).

The theorem states that $\hat{\alpha}-\alpha_0$ is asymptotically normally distributed, but its limiting distribution is centered away from zero by an $O\left(h^2 + \frac{\log \bar{N} }{\bar{n}^2h^{2d_{\zeta}}}\right)$ bias term. The first component $h^2$ arises from the bias of the imputed links and is typically the dominant term. By contrast, the second component is asymptotically negligible relative to $h^2$ in most empirically relevant settings\footnote{ To see this, note that the bandwidth minimizing the mean squared error of the imputed links is of order $h^*\sim n^{\frac{-1}{4 + d_{\zeta}}}$. For the second bias component to dominate the first, that is, $h^2 = o\left(\frac{1}{n^2h^{2d_{\zeta}}}\right)$, we would require $h = o\left(n^{-\frac{1}{1 + d_{\zeta}}}\right)$. This bandwidth shrinks much faster than the MSE-optimal bandwidth $h^*$ and therefore requires substantial undersmoothing, which is uncommon in practice. }. The asymptotic variance consists of two components: $\Omega/M$ corresponds to the asymptotic variance that would arise if the network $A$ were perfectly observed, while $V/M$ captures the additional variance introduced by nonparametric imputation of the adjacency matrix. While the bias does not necessarily shrink as $M \to \infty$, the variance component induced by network imputation decreases at rate $1/M$, this highlights an important difference between the asymptotic behavior of the downstream estimator and that of the link-level imputation error.

This result has important implications for downstream inference. Since the bias term $B$ dominates $V^{1/2}/\sqrt{M}$ in most empirically relevant settings, valid inference for $\alpha_0$ requires choosing the bandwidth in the imputation step so that $B$ is asymptotically negligible relative to the sampling variation, $\Omega^{1/2}/\sqrt{M}$. Consequently, bandwidth selection for valid downstream inference should prioritize reducing the bias term rather than balancing the bias and variance of the link-level imputation error, which requires undersmoothing in the first-stage imputation step.

In practice, we recommend that applied researchers start from the optimal bandwidth $h$ selected by cross-validation in the imputation step, then gradually decrease $h$ and examine the corresponding estimates $\hat{\alpha}$. Once the estimates become unstable as $h$ decreases, we should stop undersmoothing at that point.

Simulation Study

In this section, we examine the finite-sample properties of the proposed method using a series of numerical experiments. We compare its performance with several alternative imputation approaches, including the low-rank imputation methods proposed by bai2021matrix and li2023link, and the local PCA method developed by feng2023optimal. We also use numerical experiments to illustrate (i) the importance of incorporating observed covariates for improving prediction accuracy, and (ii) the influence of sample splitting on finite-sample performance. We begin with introducing alternative method.

Alternative Methods

\paragraph{Covariate-only method} This method is a dyadic nonparametric regression-based imputation approach that relies solely on the observed covariates $X$. When unobserved heterogeneity is present, the prediction error of this approach does not diminish as the sample size increases. We compare its numerical performance with that of our proposed method to highlight the importance of accounting for unobserved heterogeneity. We refer to this method as X in the following text.

\paragraph{Low-rank imputation} When the observed matrix follows a strong-factor structure with an additive noise component and the missing entries exhibit a block-missing pattern\footnote{ bai2021matrix refer to this pattern as a tall-wide structure in the context of panel data, while li2023link study an analogous setting in network data and refer to it as egocentric sampling. } (see (ref)), bai2021matrix and li2023link propose using low-rank estimation to impute the missing data. The key idea is to first recover the latent factors using the observed submatrices $A_{\mathcal{S}\mathcal{S}}$, $A_{\mathcal{S}\mathcal{S}^c}$, and $A_{\mathcal{S}^c\mathcal{S}}$, and then reconstruct the missing bottom-right submatrix $A_{\mathcal{S}^c\mathcal{S}^c}$ through the product of the estimated factors.

When the graphon $f$ is nonlinear as in Example (ref), the graphon matrix is typically high-rank and thus violates the low-rank assumption required by bai2021matrix and li2023link. Although their methods are not theoretically valid in such settings, we include them in our numerical experiments for comparison, given their simplicity and computational efficiency, to assess how our proposed method performs in practice relative to these benchmarks. It is worth noting that both bai2021matrix and li2023link assume that the true number of factors is known a priori. In practice, this information is rarely available. Hence, in our numerical experiments, we select the number of factors via cross-validation.

\paragraph{Local PCA} The local PCA method proposed by feng2023optimal is suitable for high-rank graphon. Specifically, it first finds the neighbors of each node $i$ based on the pseudo-distance using $k$-nearest neighbors (kNN), and then performs low-rank imputation on the submatrix formed by node $i$'s neighbors to impute the missing entries. As discussed earlier, when the underlying graphon is smoother than twice continuously differentiable, the local PCA method can, in theory, achieve a higher-order approximation error. We compare our method with local PCA to assess whether this theoretical advantage persists in finite samples.

A practical concern is that the local PCA method requires prior knowledge of the latent dimension $d_{\zeta}$, since this determines the rank when performing PCA on each submatrix. However, $d_{\zeta}$ is typically unknown and rarely supported by prior knowledge. While feng2023optimal propose an ad hoc procedure in their numerical simulations to address this issue, we instead use cross-validation to jointly select both the number of neighbors $k$ in the kNN step and the rank used in PCA.

It is worth noting that local PCA does not incorporate observed covariates\footnote{ feng2023optimal briefly discuss how to include covariates in a linear form in their appendix. }. In our numerical experiments, we therefore consider two implementations of their method: (1) LPCA, which ignores covariates $X$ and applies local PCA directly to the adjacency matrix, and (2) X-LPCA, which first removes the variation explained by $X$ and then applies local PCA to the residual matrix.

\paragraph{Local TWFE without covariates} We also consider an alternative method that directly applies our local TWFE approach to the adjacency matrix without extracting the variation explained by $X$. As discussed earlier, because the information about $\zeta$ recovered from observed data is noisy, it is unwise to ignore the observed covariates and absorb them into the latent factors. We include the numerical simulation results of this method to evaluate the importance of incorporating covariate information. We refer to this method as LTWFE and to our proposed method as X-LTWFE in the following text.

\paragraph{Sample splitting} While our asymptotic analysis relies on sample splitting for theoretical tractability, this requirement is primarily technical. In practice, sample splitting reduces the effective sample size, and since network sizes are relatively small in most economics applications, its finite-sample performance can be worse than that of the full-sample implementation. Accordingly, we implement all alternative methods without sample splitting.

For our proposed method, we consider both versions, with and without sample splitting, to examine how this choice affects finite-sample performance. We refer to our method with sample splitting as X-LTWFE-SP in the following text.

\paragraph{Summary} We summarize all imputation methods used in the numerical experiments below:

Simulation for Imputation Accuracy

The network in the simulation is generated as

The dimensions of the observable and latent characteristics are $d_{X} = d_{\xi} = 2$, and $\{(X_i, \xi_i)\mid i=1,\ldots, N\}$ is i.i.d. across individuals. Covariates and latent factors are generated according to

and

The idiosyncratic shocks $\{U_{ij}\}_{1 \leq i < j \leq N }$ are independently drawn from a standard logistic distribution. The coefficient vector $\beta$ captures homophily effects on observables and takes values $\beta\in \left\{(-0.5, -0.5), (-2, -2)\right\}$. We vary the magnitude of $\beta$ to control the sparsity of the network: for example, when $\beta = (-0.5, -0.5)$, the resulting network tends to be denser, with a larger number of links being formed.

The network size is $N = 200$, and in the simulations we consider sampling rates $\varphi := n/N \in \{0.2, 0.3, 0.4, 0.5, 0.6, 0.8\}$ to compare the performance of the proposed and alternative methods under different sampling scenarios. We conduct $S = 1,000$ Monte Carlo replications in the numerical experiment, and for each replication $s$, we construct the imputed network $\hat A_s$ and compute the mean squared error (MSE) between $\hat A_s$ and the underlying probability $P_{s}$ on the missing subnetwork, $ \frac{1}{(N-n)(N-n-1) } \sum_{i, j\in \mathcal{S}^c} \left(\hat{A}_{s, ij} - P_{s, ij} \right)^2$. We report the root mean-squared error (RMSE)

We evaluate the performance of our imputation method and alternative approaches in both dense networks (e.g., $\beta = (-0.5, -0.5)$, for which the average probability of forming a link is about $34\%$) and sparse networks (e.g., $\beta = (-2, -2)$, for which the average probability of forming a link is about $10\%$). The simulation results are reported in Table (ref), with RMSEs presented in units of $0.01$ for ease of reading.

The first column (X) reports the RMSE of covariate-only method, which imputes missing links using $X$ alone. Because this method ignores unobserved heterogeneity, its imputations are systematically biased, and the error does not vanish as the network size or sampling rate increases. Consistent with the theoretical prediction, the numerical results show that, as the sampling rate $\varphi$ increases, the RMSE of the covariate-only method exhibits very limited improvement (from $0.210$ to $0.206$ in the dense scenario, and from $0.138$ to $0.135$ in the sparse scenario).

The second column (LR) reports the RMSE of the low-rank imputation, and the third column (LPCA) reports the RMSE of the local PCA. Since the underlying graphon in (ref) is high-rank, local PCA should, in theory, outperform the low-rank approach. However, the numerical results indicate that the RMSEs of the two methods are nearly the same across all scenarios. We conjecture that this pattern is driven by two reasons: (1) Although the graphon is not exactly low-rank, it is approximately low-rank, so the low-rank approach still provides a good approximation to the graphon matrix, especially when performance is measured using Frobenius norm. (2) Local PCA relies only on the subnetwork formed by neighbors, which leads to substantial sample loss and limits its performance in finite samples. In line with our theory, the local PCA estimator outperforms the covariate-only (X) method in most cases, and its prediction error decreases as the sampling rate increases.

The fourth column (LTWFE) reports the RMSE of the local two-way fixed0effects imputation that ignores covariates. In line with the theory, its prediction error decreases as $n$ increases. Table (ref) shows that at low sampling rates, the LTWFE method performs better than local PCA, whereas once the sampling rate becomes large (e.g., when $\varphi$ exceeds $40\%$), local PCA outperforms LTWFE. This suggests that although local PCA has better asymptotic properties, in the sense that it achieves higher-order approximation error, it requires a larger sample size to realize this advantage, whereas LTWFE is more robust in small samples. In addition, the LTWFE method is more robust to network sparsity: as shown in Table (ref), its performance is nearly unaffected by changes in the sparsity of the adjacency matrix.

The fifth column (X-LPCA) and the sixth column (X-LTWFE) report the performance of LPCA and LTWFE when incorporating covariate information. These methods first extract the variation of the missing links that can be explained by $X$, and then apply LPCA and LTWFE to the residual matrix, respectively. Compared with LPCA and LTWFE, incorporating covariate information leads to substantial improvements in prediction accuracy, particularly for X-LTWFE. The gains are especially large when the sampling rate is low or when the network is sparse. X-LPCA exhibits a similar pattern. In addition, Table (ref) shows that the proposed method outperforms all alternatives in our simulations, and its performance remains robust even under sparse networks and low sampling rates, which are common in empirical applications.

The last column (X-LTWFE-SP) reports the performance of our proposed estimator with sample splitting. As discussed earlier, although sample-splitting facilitates asymptotic analysis, it reduces the effective sample size, and this issue can be severe in small-sample settings. In addition, it is reasonable to expect (though we do not formally prove it) that the correlation between the pseudo-distance and the error terms diminishes, making sample splitting unnecessary in theory. Table (ref) shows that performance of proposed imputation under sample splitting is always worse than the performance without splitting the sample, but the difference between the two methods becomes smaller as the sample size increases.

In summary, the simulations show that our proposed method outperforms the alternative methods and is more robust across different levels of network sparsity and sampling rates. In addition, for finite samples, we recommend using the version without sample splitting in order to achieve higher prediction accuracy. Further numerical results are provided in Appendix (ref).

Simulation for GMM estimators

In this subsection, we evaluate the performance of GMM estimators in downstream analyses that rely on the imputed networks. We compare the performance of our method with the alternative approaches introduced previously. We focus on two widely used classes of empirical exercises in economics that are based on network data. \paragraph{Regression on network statistics} Suppose for each network $m=1, \ldots, M$, the outcome $Y_{mi}$ is generated according to:

where $\phi_{mi} : [0,1]^{N_m \times N_m} \to \mathbb{R}$ computes a network statistic for individual $i$ in network $m$, $\{u_m\}_{1\leq m\leq M}$ are i.i.d. network-level random effects independent of $A_m$, and $e_{mi}$ is i.i.d. idiosyncratic error. We consider two specifically, normalized average degree and eigenvector centrality, as defined in Example (ref). The networks $\{A_m\}_{1\leq m\leq M}$ according to (ref) with $\beta = (-0.5, -0.5)$, which yields relatively dense networks with a larger number of links.

We set the number of networks to $M = 40$ and fix network size at $N_m = 200$ for all $m$. The simulations consider sampling rates $\varphi := n/N \in \{0.2, 0.3, 0.4, 0.5\}$ to compare the performance of the proposed method and alternative methods under different levels of sampling. For reference, we also report the performance of regressions using the complete networks (CD). We set $\alpha_0 = 0$ and $\alpha_1 = 0.5$.

We conduct $S = 1,000$ Monte Carlo replications. For each replication $s$, we generate $M$ independent network, impute each network, and estimate the linear regression. The bias and standard deviation of $\hat{\alpha}_1$ are reported in Table (ref).

Table (ref) displays patterns similar to those in Table (ref). As the sampling rate increases, both the bias and the standard deviation of the estimator decline in most cases. This improvement is not only due to more accurate imputation at higher sampling rates (Theorem (ref)), but also because the aggregate statistics themselves are computed more precisely when fewer observations are missing. The latter explains why the bias and standard deviation in the X column also decrease. In most cases, the X-LTWFE estimator delivers the best performance and remains stable even when the sampling rate is small. The only exception arises when eigenvector centrality is used and the sampling rate is high. In this case, LTWFE yields an extremely small bias—almost identical to that obtained using the complete data. This pattern appears to be driven by features of our numerical design.

It is also worth noting that the bias-std ratio of the linear regression estimator is not close to zero in these simulations. This arises from the use of nonparametric imputation and is particularly evident when eigenvector centrality is employed. Even if researchers choose to rely on a parametric model for imputation, the robustness of our estimator implies that our method remains valuable as a robustness check in applied work.

\paragraph{Linear-in-mean peer-effects model} Suppose for each network $m=1, \ldots, M$, the outcome $Y_{mi}$ is generated according to:

where $\{u_m\}_{1\leq m\leq M}$ are i.i.d. network-level random effects independent of $A_m$, and $e_{mi}$ is i.i.d. idiosyncratic error. The networks $\{A_m\}_{1\leq m\leq M}$ according to (ref) with $\beta = (-0.5, -0.5)$, which yields relatively dense networks with a larger number of links.

We set the number of networks to $M = 40$ and fix network size at $N_m = 200$ for all $m$. The simulations consider sampling rates $\varphi := n/N \in \{0.2, 0.3, 0.4, 0.5\}$ to compare the performance of the proposed method and alternative methods under different levels of sampling. For reference, we also report the performance of regressions using the complete networks (CD). We set $\alpha_C = 0$, $\alpha_{\bar{Y}} = 0.5$, $\alpha_W = \alpha_{\bar{W}} = (1, 1)$. We conduct $S = 1,000$ Monte Carlo replications. For each replication $s$, we generate $M$ independent network, impute each network, and estimate the GMM estimator using (ref) with $\Sigma_M = \boldsymbol{I}_{q}$. The bias and standard deviation of estimate of endogenous effect $\hat{\alpha}_{\bar{Y}}$ are reported in Table (ref).

Table (ref) reports the bias and standard deviation of different estimators of $\alpha_{\bar{Y}}$ and $\alpha_{W,1}$ obtained using X, LTWFE, X-PCA, X-LTWFE, and sample-splitting X-LTWFE. The full set of simulation results is presented in Table (ref) and Table (ref) in Appendix (ref). Similar to the results in Table (ref), as the sampling rate increases, both the bias and the standard deviation of the estimator decline. In most cases, the estimator based on our proposed method (X-LTWFE) has strong performance, especially when the sampling rate is small. However, we also observe that, when the sampling rate exceeds $30\%$, X-LPCA delivers strong performance and attains a smaller bias in estimating the endogenous effect $\alpha_{\bar{Y}}$ than our proposed estimator. In addition, although LTWFE delivers a smaller bias in estimating the endogenous effect than X-LTWFE when $\varphi \geq 40\%$, its performance deteriorates substantially when the sampling rate becomes small.

The numerical patterns align with the implications of Theorem (ref): the bias-to-std ratio of the estimators does not vanish to $0$. This phenomenon is inevitable when the graphon is a nonparametric function of latent factors and poses challenges for inference. Nevertheless, by Theorem (ref), our proposed estimator is still valuable as a robustness benchmark in empirical applications.

Empirical Application

This section presents an empirical application using the network data of 43 villages in rural India from banerjee2013diffusion. For each village, the researchers informed certain households (injection nodes or leaders) about the information of microfinance program and asked them to encourage other villagers to join. The researchers want study how households' microfinance take-up decisions are influenced by network characteristics. The data consist of egocentrically sampled social network information and demographic variables (for both sampled and unsampled households). There are $M = 43$ villages, with an average of approximately $N \approx 200$ households per village, and a sampling rate of about $\varphi = 45\%$.

We need to impute the missing links in the adjacency matrices due to egocentric sampling. Since the dataset includes rich demographic information with $d_X > 10$ and the number of households within each village is relatively small, fully nonparametric regression becomes infeasible. Thus, we replace the first-step nonparametric regression with a linear projection on $\omega(X_i, X_j): = ((X_{i1} - X_{j1})^2, \ldots, (X_{id_{X}} - X_{jd_{X}})^2)$ (another option is to select a subset of important variables and then conduct a nonparametric regression, but this would rely heavily on researchers' discretion). We then apply the local two-way fixed-effects regression to the residuals from the projection.

We focus on the household-level regressions in banerjee2013diffusion. The outcome variable is each household's microfinance take-up decision ($Y_i = 1$ if the household decides to take part in the program). We consider three specifications: (i) how households' normalized average degrees (defined as in Example (ref)) affect their take-up decisions; (ii) how households' eigenvector centralities (defined as in Example (ref)) affect their take-up decisions; and (iii) how peers' decisions affect individual take-up through a linear-in-means peer-effects model.

Table (ref) reports estimates of household-level regressions in which the dependent variable is a household's microfinance take-up decision and the regressors are household-level network characteristics. We provide regression results for both sampled (incomplete) networks and imputed network using our proposed method. Panel A shows that, when using the imputed network, the household's take-up probability increases by $7.2\%$ percentage points when its normalized average degree increases by $10\%$ (equivalently, when the household has about $20$ more connections within the village). The estimated effect is weaker when using sampled networks, with the coefficient $\alpha_1$ decreasing from $0.72$ to $0.61$. Panel B presents results using eigenvector centrality. Both specifications yield similar estimates ($0.026$ using the sampled network and $0.028$ using the imputed network), indicating that households with higher eigenvector centrality are more likely to take up microfinance.

Table (ref) reports estimates from a household-level linear-in-means peer-effects model, where the dependent variable is the household's microfinance take-up decision, and the regressor $W$ is a dummy variable indicating whether the household contains an injection node. We present regression results based on both the sampled (incomplete) network and the imputed network constructed using our proposed method.

The coefficient $\alpha_{\bar{Y}}$ captures the endogenous peer effect, measuring the influence of neighbors' participation decisions. Using the imputed network, the estimate of $\alpha_{\bar{Y}}$ is $0.09$, implying that a household's take-up probability increases by approximately $0.9$ percentage points when the fraction of its friends who join the program rises by $10$ percentage points (corresponding to roughly two additional friends joining the program, on average). In contrast, the estimated endogenous effect is $-0.02$ negative when using sampled networks. In addition, both regressions give similar estimates for the direct effect, suggesting that the household's take-up probability increases by $6$ percentage points when it contains an injection individual. Finally, for the contextual peer effect $\alpha_{\bar{W}}$, both sampled and imputed networks produce negative estimates, with the magnitude being larger when using the imputed network.

Conclusion

Sampled network data are common in empirical research because collecting full network information is costly, but using sampled networks can lead to biased estimates. We propose a flexible imputation method for sampled networks and show that downstream empirical analysis using the imputed network yields consistent parameter estimates.

Our method separately imputes the part of each missing link that is explained by covariates and the remaining variation not captured by them: the former via projection onto observed covariates, and the latter using local two-way fixed effects. The imputation method avoids parametric assumptions, does not rely on low-rank restrictions, and flexibly incorporates observed covariates and unobserved heterogeneity. We establish entrywise convergence rates for the imputed matrix and prove consistency of GMM estimators based on the imputed network. We further derive the convergence rate for the corresponding estimator in the linear-in-mean peer-effects model. Simulations show strong performance of our method both in terms of imputation accuracy and in downstream empirical analysis.

The robustness of our imputation comes at the cost of making downstream inference challenging. In the paper, we discuss how undersmoothing can help improve the downstream estimator, and other approaches, such as bias correction may be potentially applicable. Developing inference procedures for downstream estimators is a potential topic for future research.