Graph Neural Networks for Causal Inference Under Network Confounding

Graph Neural Networks for Causal Inference Under Network Confounding

\onehalfspacing

\addcontentsline{toc}{part}{Main Paper}

Introduction

Treatment assignment is said to be unconfounded if it is as good as random within subpopulations of observationally equivalent units. When the stable unit treatment value assumption (SUTVA) is plausible, units with identical covariates are naturally considered observationally equivalent. However, when units are connected through a network, they may differ in other observed dimensions that may confound causal inference under network interference. These dimensions may include, for example, the number of type-$x$ neighbors, the number of type-$x$ neighbors with $m$ neighbors of type $y$, and so on through higher-order neighbors.

Existing work typically adjusts for confounding only on a small subset of these dimensions such as the average covariates of neighbors. However, it may be difficult to justify a particular choice in practice. Neighbor covariates can influence selection into treatment in more complex ways not adequately captured by the mean or other statistic of convenience.

This paper considers a more general formulation of unconfoundedness. Our objective is estimation of and inference on reduced-form measures of treatment and spillover effects using data from a single large network. We model selection as the outcome of a simultaneous equations model, which allows for peer effects. As a result, treatments are functions of the entire network $\bm{A}$ and all unit covariates $\bm{X}$, and it is generally not possible to summarize confounding by a simple low-dimensional function of these objects. Our unconfoundedness condition instead considers units observationally equivalent if they occupy identical positions in the network, meaning that they match on all observed neighborhood and higher-order neighborhood characteristics.

Methods that rule out interference in selection may result in biased estimates. For example, consider the causal effect of vaccine adoption on illness. With peer effects in adoption, vaccinated individuals tend to have more vaccinated direct and indirect social contacts, and a simple comparison of adopters and non-adopters may overstate vaccine efficacy, even after controlling for neighbors' covariates.

Because peer effects in selection induce high-dimensional network confounding, they are challenging to accommodate in a nonparametric setting. Most of the literature on interference, including recent work using instrumental variables for identification, rules out selection peer effects ditraglia2023identifying,hoshino2023causal,kang2016peer. Papers that allow for them typically rely on semiparametric, game-theoretic models of selection, which substantially reduce the dimensionality of the problem but may be subject to model misspecification hoshino2023treatment,jackson2020adjusting,kim2020analysis,lin2021selection.\footnote{Two exceptions are balat2023multiple, who study partial identification under strategic complementarities, and imai2021causal, who study two-stage randomized trials with noncompliance. Both consider nonparametric selection models but require data consisting of many independent networks.}

To account for network confounding, an initial idea might be to apply double machine learning using the lasso for high-dimensional estimation chernozhukov2018double. Under SUTVA, implementation of the lasso requires the specification of a basis $\{P_k(X_i)\}_{k=1}^d$ for unit-level covariates $X_i$. However, in our setting, a unit $i$'s “covariates” correspond to its network position $(i,\bm{X},\bm{A})$ and it is unclear how to choose a basis $\{P_k(i,\bm{X},\bm{A})\}_{k=1}^d$ for such an object. Nothing in the standard toolbox for high-dimensional estimation suggests that it is possible to nonparametrically adjust for network confounding of this sort, which is presumably why it has been avoided in the literature.

Contributions

Our primary contribution is to propose a nonparametric model with endogenous peer effects and make the case that estimation of common reduced-form estimands remains feasible despite high-dimensional network confounding. We utilize the model of approximate neighborhood interference (ANI) proposed by leung2022causal, which posits that interference decays with network path distance. Leung shows that ANI allows for endogenous peer effects but focuses on randomized assignment. In observational settings, there may also be peer effects in selection. We therefore relax his assumption of independent treatment assignment to allow for ANI in treatment selection. Since both stages are simultaneous equations models, this creates significant new complications in the form of high-dimensional network confounding.

We observe that graph neural networks (GNNs) can provide a flexible nonparametric basis for functions of network position. We propose to use doubly robust estimators with first-stage nuisance functions approximated by GNNs. Whereas a variety of conventional machine learners can be employed for the first stage under SUTVA, none are well suited to our setting, and the observation that GNNs can fill this gap is novel.

Since the first-stage nuisance functions are nonparametric, high-dimensional functions of network position, some form of low-dimensional structure is required for estimation to be feasible. Additionally, GNNs have been found empirically to perform best with shallow architectures, which correspond to relatively low-dimensional parameterizations, in contrast to deep architectures popular with convolutional neural networks alon2021bottleneck,li2018deeper. Our key insight is to draw a novel connection between ANI and approximate sparsity conditions in the lasso literature. We argue that this connection justifies the use of shallow GNN architectures in our setting.

To understand the idea, let $\mathcal{N}(i,K)$ denote $i$'s {\em $K$-neighborhood}, the set of units whose path distance from $i$ is at most $K$, and $(\bm{X}_{\mathcal{N}(i,K)}, \bm{A}_{\mathcal{N}(i,K)})$ denote the restriction of $(\bm{X},\bm{A})$ to $\mathcal{N}(i,K)$. The key parameter of a GNN is its {\em depth} or number of layers $L$, which determines the {\em receptive field} $(\bm{X}_{\mathcal{N}(i,L)}, \bm{A}_{\mathcal{N}(i,L)})$ used to predict $i$'s outcome. For example, a one-layer GNN only uses $i$'s 1-neighborhood $(\bm{X}_{\mathcal{N}(i,1)}, \bm{A}_{\mathcal{N}(i,1)})$ to predict its outcome, rather than the entirety of $(\bm{X},\bm{A})$. Accordingly, the choice of $L$ depends on prior information about the function being estimated. Under high-dimensional network confounding, it would seem that the propensity score depends nontrivially on the entirety of $(\bm{X},\bm{A})$, which would require a relatively large choice of $L$.

Under ANI, however, interference decays with distance, which implies that a unit $i$'s outcome and treatment are primarily determined by $(\bm{X}_{\mathcal{N}(i,L)}, \bm{A}_{\mathcal{N}(i,L)})$ for relatively small $L$. This is reminiscent of approximate sparsity, under which the regression function primarily depends on a small subset of regressors. As a result, our first-stage nuisance functions can be approximated by lower-dimensional analogs that only depend on the $L$-neighborhood, and these can be directly estimated with shallow $L$-layer GNNs. Our formal result provides primitive conditions on interference rationalizing small choices of $L$ of order $\log n$.

Existing work on network interference assumes the correlation structure is characterized by a dependency graph emmenegger2022treatment,ogburn2022causal, which requires dependence between units to be zero after a known network distance. This is analogous to $m$-dependence and incompatible with endogenous peer effects. We establish that the data in our setting is $\psi$-dependent kojevnikov2021limit, a correlation structure more similar to mixing and near-epoch dependence in that it allows correlation to decay smoothly with distance. Under high-level conditions on the GNNs, we show that the doubly robust estimator is approximately normally distributed in large networks. We propose a HAC-type variance estimator based on kojevnikov2021limit and kojevnikov2021bootstrap which we prove is non-negative and consistent. We also provide a bandwidth that adjusts for estimation error in the first-stage machine learners.

We substantiate the theory in a simulation study and empirical illustration revisiting the microfinance diffusion application of he2024measuring. We show how our estimands can capture different dimensions of diffusion complementary to their “average diffusion at the margin” measure. Our theoretical framework allows for more complex diffusion processes without requiring the econometrician to prespecify the maximum number of within-period rounds of diffusion. By including richer controls that account for network confounding, we find more attenuated diffusion effects.

Related Literature

There is a large literature on interference primarily focusing on randomized control trials athey2018exact,li2022random,toulis2013estimation. We contribute to a growing recent literature on unconfoundedness, much of which operates in a partial interference setting where units are partitioned into disjoint groups with no interference across groups liu2019doubly,qu2022efficient.

veitch2019using consider network interference and propose to use as controls “node embeddings,” which are learned functions of the graph. Node embeddings can be obtained from a variety of methods, but it is unclear how to justify a particular choice. GNNs learn node embeddings in an end-to-end fashion (see (ref)), and our behavioral model provides justification for their use.

auerbach2022identification studies identification conditions distinct from unconfoundedness but proposes a related strategy of “matching” on certain network statistics. He provides conditions under which pairwise differencing using unit pairs matched on a novel codegree statistic eliminates selection bias.

Prior to the GNN literature, graph kernels were the dominant method for graph learning tasks morris2021weisfeiler. These are to kernel regression as GNNs are to sieve estimation, so graph kernels require a user-specified measure of similarity between regressors, in this case, between two graphs. auerbach2023local propose a graph kernel estimator using a novel similarity measure based on graph isomorphism. Since there is no known algorithm for isomorphism testing with polynomial runtime in the network size ((ref) discusses GNNs' relationship to this problem), many graph kernel approaches amount to specifying an “embedding,” a mapping from networks to Euclidean space kriege2020survey. As noted by wu2020comprehensive, embeddings are predetermined functions of the network, whereas GNNs produce learnable embeddings.

Finally, our paper contributes to a recent econometric literature using neural networks as sieve estimators for regression functions athey2021using,farrell2021deep,kaji2020adversarial. Whereas other machine learners can theoretically be employed in these settings, our problem cannot be solved with standard methods. In this respect, our work relates to pollmann2021causal who employs convolutional neural networks to construct counterfactuals for spatial treatments.

The next section defines the model and estimators. We introduce GNNs in (ref) and characterize the asymptotic properties of our estimators in (ref). In (ref), we provide conditions under which the nuisance functions have low-dimensional structure. Section (ref) reports results from a simulation study, and (ref) presents the empirical application. Finally, (ref) concludes.

We represent an undirected network $\bm{A}$ as an $n\times n$ binary adjacency matrix with $ij$th entry $A_{ij} \in \{0,1\}$ representing a link between units $i$ and $j$. We assume no self-links, so $A_{ii}=0$. Let $\ell_{\bm{A}}(i,j)$ denote the {\em path distance} between $i,j$ in $\bm{A}$, defined as the length of the shortest path between them, if one exists, and $\infty$ otherwise. The {\em $K$-neighborhood} of a unit $i$ in $\bm{A}$ is denoted by $\mathcal{N}(i,K) = \{j\in\mathcal{N}_n\colon \ell_{\bm{A}}(i,j) \leq K\}$ and its size by $n(i,K) = \lvert\mathcal{N}(i,K)\rvert$. We refer to the elements of $\mathcal{N}(i,K)\backslash\{i\}$ for $K=1$ as $i$'s {\em neighbors} and the elements of the same set for $K>1$ as $i$'s {\em higher-order} neighbors. A unit $i$'s {\em degree} is $n(i,1)-1$, the number of neighbors.

Setup

Let $\mathcal{N}_n = \{1, \dots, n\}$ be the set of units connected through the network $\bm{A}$. Each unit $i\in\mathcal{N}_n$ is endowed with unobservables $(\varepsilon_i, \nu_i) \in \mathbb{R}^{d_\varepsilon} \times \mathbb{R}^{d_\nu}$ and observables $X_i \in \mathbb{R}^{d_x}$. The model primitives determine outcomes and treatments according to

respectively, where $\bm{X} = (X_i)_{i=1}^n$ is the matrix with $i$th row equal to $X_i'$; $\bm{Y}$, $\bm{D}$, $\bm{\varepsilon}$, and $\bm{\nu}$ are similarly defined; and $\{(g_n,h_n)\}_{n\in\mathbb{N}}$ is a sequence of function pairs such that each $g_n(\cdot)$ has codomain $\mathbb{R}$ and $h_n(\cdot)$ codomain $\{0,1\}$. The econometrician observes $(\bm{Y}, \bm{D}, \bm{X}, \bm{A})$, and our analysis treats $(\bm{A}, \bm{X}, \bm{\varepsilon}, \bm{\nu})$ as random.\footnote{The asymptotic theory in (ref) conditions on $(\bm{X},\bm{A})$ to avoid imposing additional assumptions on its dependence structure. A design-based analysis would additionally condition on $\bm{\varepsilon}$, but this would generally preclude consistent estimation of the nonparametric functions in the doubly robust estimator defined in (ref).}

The timing of the model is as follows. First, nature draws the primitives $(\bm{A},\bm{X},\bm{\varepsilon},\bm{\nu})$. Next, units select into treatment, potentially based on other units' decisions, and $h_n(\cdot)$ is the reduced-form outcome of that process. Finally, $g_n(\cdot)$ is the reduced form of the subsequent process that generates outcomes. Because $g_n(\cdot)$ and $h_n(\cdot)$ may depend on the primitives of all units, the setup allows $\bm{Y}$ and $\bm{D}$ to be outcomes of simultaneous equations models with endogenous peer effects, as shown in the next examples.

Given specification (ref), we define potential outcomes as

Confounding may arise because $Y_i(\bm{d})$ is potentially correlated with $D_i$ due to high-dimensional observables $(\bm{X},\bm{A})$ and the unobservables $\bm{\varepsilon}$ and $\bm{\nu}$ may be correlated. We restrict the second source of confounding.

As discussed below, unconfoundedness conditions used in the existing literature additionally limit the first source of confounding to known summary statistics of $(\bm{X},\bm{A})$. Ours is more analogous to standard formulations of unconfoundedness under SUTVA since we do not impose an index restriction on observed confounding.

Because the econometrician only observes a single network, a large-sample theory requires restrictions on interference to obtain some form of weak dependence. We next specify a nonparametric model of decaying interference that accommodates the previous examples. For any $S \subseteq \mathcal{N}_n$, let $\bm{D}_S = (D_i)_{i\in S}$, and similarly define $\bm{X}_S$ and other such submatrices. Let $\bm{A}_S = (A_{ij})_{i,j \in S}$ denote the subnetwork of $\bm{A}$ on $S$, formally the submatrix of $\bm{A}$ restricted to $S$. Recall that $\mathcal{N}(i,s)$ is the $s$-neighborhood of $i$ in $\bm{A}$.

This is analogous to the model of approximate neighborhood interference proposed by leung2022causal but imposed on both the outcome and selection models. Whereas $g_n(i,\dots)$ is unit $i$'s realized outcome, $g_{n(i,s)}(i,\dots)$ is its outcome under a counterfactual “$s$-neighborhood model.” In the latter case, we fix all model primitives and treatments at their realized values, drop units outside of $\mathcal{N}(i,s)$ from the model, and direct the remaining units to interact according to the process $g_{n(i,s)}(\cdot)$ to produce counterfactual $s$-neighborhood outcomes.\footnote{This formulation of ANI is related to an estimation strategy proposed by xu2018social for binary games on networks with incomplete information. His idea is to approximate an agent $i$'s strategy in the $n$-agent game with its strategy in the counterfactual game restricted to $i$'s $s$-neighborhood.} The error from approximating the observed outcome with the $s$-neighborhood counterfactual is bounded by $\gamma_n(s)$, which decays with the neighborhood radius $s$. This formalizes the idea that $Y_i$ is primarily determined by units relatively proximate to $i$, so that the further a unit is from $i$, the less it influences $i$'s outcome. The second equation imposes the analogous requirement on $D_i$. It is trivially satisfied if $D_i$ is i.i.d., in which case $h_n(i, \bm{X}, \bm{A}, \bm{\nu})$ is only a function of $\nu_i$ for i.i.d.\ $\{\nu_i\}_{i=1}^n$.

Previous Approaches

The first contribution of this paper is to propose a potential outcomes model and formulation of unconfoundedness with richer microfoundations relative to the literature. The second contribution is to show that inference remains possible despite the high-dimensional network confounding unique to our setting. We next elaborate on the first contribution by comparing our model to the literature. The second contribution will be the focus of the remainder of the paper.

The standard SUTVA model and unconfoundedness condition correspond to

To generalize this setup to allow for network interference, the literature proceeds as follows. Define

where $f_n(\cdot)$ and $q_n(\cdot)$ are known vector-valued functions. The {\em effective treatment} manski2013identification or {\em exposure mapping} aronow_estimating_2017 $T_i$ is a low-dimensional function of the treatment assignment vector. The {\em network controls} $W_i$ are low-dimensional functions of the covariates. The literature commonly employs the {\em neighborhood interference} model and unconfoundedness condition

which is a direct generalization of (ref) emmenegger2022treatment,forastiere2021identification,ogburn2022causal. Here $T_i$ entirely summarizes interference while $W_i$ summarizes confounding.

Common examples of $T_i$ and $W_i$ are

Versions of these are employed in the numerical illustrations of the previous three references. The specification of $T_i$ implies that $Y_i$ depends on $\bm{D}$ only through two statistics: own treatment and the number of treated neighbors. Under model (ref), variation in the first component identifies a direct treatment effect and variation in the second a spillover effect. Like most exposure mappings used in the literature, this only depends on $\bm{D}_{\mathcal{N}(i,1)}$, so the outcome model (ref) implies no interference beyond the 1-neighborhood. Likewise, the choice of $W_i$ implies no confounding beyond 1-neighborhood covariates. More generally, one could restrict the outcome model to depend only on the $K$-neighborhood treatments $\bm{D}_{\mathcal{N}(i,K)}$ for some fixed threshold $K$. Unlike (ref), this does not allow for endogenous peer effects.

Whereas leung2022causal focuses on randomized experiments, we study observational data on economic agents who choose to select into treatment. It then becomes important to specify a model of selection rationalizing the econometrician's choice of controls. sanchez2022spillovers is the first to provide such a model. Under neighborhood interference (ref) and an exposure mapping similar to (ref), he shows that it is sufficient to set $W_i=X_i$, that is, to solely control for own covariates. Since much of the literature utilizes controls such as (ref), this raises the question of what model of selection justifies their use or more broadly the use of “network controls” that depend more generally on $(\bm{X},\bm{A})$.

Our model (ref) provides an answer. The presence of complex interference in both the outcome and treatment stages induces selection on $(\bm{X},\bm{A})$, so that it is generally insufficient to control only for a simple summary statistic such as (ref). Our outcome model is considerably more general than (ref) because we do not require existence of a low-dimensional function $T_i$ of $(\bm{D},\bm{A})$ summarizing interference. Our unconfoundedness condition ((ref)) is likewise considerably weaker than (ref) because we do not require existence of a low-dimensional function $W_i$ of $(\bm{X},\bm{A})$ summarizing confounding.

Estimand

Following much of the literature, we focus on estimands defined by exposure mappings leung2024causal,savje2021causal, though the core idea of accounting for high-dimensional network confounding using GNNs may potentially be applied to other estimands. Recall from the previous subsection the definition of the exposure mapping $T_i = f_n(i, \bm{D}, \bm{A})$, where $\{f_n\}_{n\in\mathbb{N}}$ is a sequence of functions with common codomain $\mathcal{T}$, a discrete subset of $\mathbb{R}^{d_t}$. Let $\mathcal{M}_n \subseteq \mathcal{N}_n$ be a subset of the units and $m_n = \lvert\mathcal{M}_n\rvert$. We study the estimand

This compares average outcomes of units under two different values of the exposure mapping while adjusting for high-dimensional network confounders. The comparison is restricted to a subpopulation $\mathcal{M}_n$, the choice of which can be important for ensuring overlap, as discussed in the following examples. Depending on the choice of $f_n(\cdot)$, $t$, and $t'$, $\tau(t,t')$ may be intended to measure an average treatment or spillover effect, as illustrated in the examples below.

Our large-sample results pertain to the following subpopulations and exposure mappings, which include the previous examples.

In (ref) for $t=(0,1)$, this holds for $d=0$, $\bm{\Delta} = (0,\infty)$, and $\bm{\Gamma}$ given in the example. In (ref) with $t=(0,2)$, this holds for $d=0$, $\bm{\Delta}=[1.5,2.5]$, and $\bm{\Gamma}=[2.5,3.5]$.

We restrict to this class of mappings for two reasons. First, it includes the most widely used examples in the literature, which are those presented above. Second, $T_i$ is a complex and discontinuous function of treatments, which are in turn complex functions of the structural primitives. We require additional structure on $T_i$ to characterize the dependence structure of the data and apply a CLT.

leung2024causal provides conditions under which $\tau(t,t')$ has a causal interpretation. The following result is a direct consequence of his Theorem A.1. Let $p_{i,t}(\cdot \mid \bm{x},\bm{a})$ denote the conditional probability mass function of $\bm{D} \mid T_i=t, \bm{X}=\bm{x}, \bm{A}=\bm{a}$.

The result represents $\tau(t,t')$ as a convex average of unit-level effects of the form $Y_i(\bm{d}_{i,t}, \bm{d}_{-\mathcal{N}(i,1)}) - Y_i(\bm{d}_{i,t'}, \bm{d}_{-\mathcal{N}(i,1)})$ such that the exposure mapping evaluates to $t$ under the assignment vector $(\bm{d}_{i,t}, \bm{d}_{-\mathcal{N}(i,1)})$ and likewise for $t'$. In (ref) with $t=(0,2)$, $\bm{d}_{i,t}$ is the vector with dimension $n(i,1)$ such that $i$ is untreated and has exactly two treated neighbors. Furthermore, $\bm{d}_{i,t},\bm{d}_{i,t'}$ are partially ordered, which means that increasing the exposure mapping from $t$ to $t'$ pushes more neighborhood units into treatment.

Treatments are conditionally independent if treatment adoption follows a nonparametric game of incomplete information $D_i = h_n(i, \bm{X}, \bm{A}, \nu_i)$ as in (ref) and private information $\nu_i$ is independently distributed across units conditional on $(\bm{X},\bm{A})$. Structural analyses commonly assume private information is i.i.d.\ and independent of public information lin2021selection,xu2018social.

Estimator

Define the nuisance functions

which we refer to respectively as the (generalized) propensity score imbens2000role and outcome regression. Let $\hat{p}_t(i,\bm{X},\bm{A})$ and $\hat{\mu}_t(i,\bm{x},\bm{a})$ denote their respective GNN estimators, which will be defined in (ref). We use a standard doubly robust estimator for multi-valued treatments

where

To define the variance estimator, we need some notation. Define the bandwidth

where $\lceil\cdot\rceil$ rounds up to the nearest integer, $\delta(\bm{A}) = n^{-1} \sum_{i=1}^n (n(i,1)-1)$ is the average degree, and $\mathcal{L}(\bm{A})$ is the average path length.\footnote{We assume $\delta(\bm{A})>1$, as is typical in practice. By the average path length, we mean the average over all unit pairs in the largest component of $\bm{A}$. A component is a connected subnetwork such that all units in the subnetwork have infinite path distance to non-members of the subnetwork.} Define $\bm{K}^U$ and $\bm{K}^{PD}$ as $m_n\times m_n$ matrices with respective $ij$th entries $\mathbf{1}\{\ell_{\bm{A}}(i,j)\leq b_n\}$ and

where the rows and columns of both matrices range over $i,j \in \mathcal{M}_n$. We propose the estimator

This takes the larger of two HAC estimators using different kernels. Since $\bm{K}^{PD}$ is positive semidefinite, $\hat\sigma^2$ is always non-negative.\footnote{Let $\bm{G}$ be the matrix obtained from $\bm{K}^U$ by dividing each row of the latter by the square root of its row sum and replacing $b_n$ with $b_n/2$. Then $\bm{K}^{PD} = \bm{G}'\bm{G}$, which is positive semidefinite.}

Graph Neural Networks

Consider the problem of estimating an unknown scalar function of network position $F^*(i, \bm{X}, \bm{A})$. A GNN estimator of $F^*(\cdot)$ is a parameterized function that maps $(\bm{X},\bm{A})$ to a vector of unit-level predictions $(\hat{F}(i, \bm{X}, \bm{A}))_{i=1}^n \in \mathbb{R}^n$. In (ref) we define the standard GNN architecture. In (ref), we present GNN estimators for the nuisance functions. In (ref), we discuss the significance of a nonparametric shape restriction, permutation invariance, that GNNs impose.

Architecture

The standard GNN architecture consists of nested, parameterized, vector-valued functions called {\em neurons} that are arranged in $L$ {\em hidden layers} with $n$ neurons per layer. Let $h_i^{(l)}$ denote the $i$th neuron in hidden layer $l$, which is typically a vector of dimension or {\em width} $H$ that is common across $i$ and $l$. This is often interpreted as unit $i$'s {\em node embedding}, a Euclidean representation of its network position. As we progress to higher-order layers, say $h_i^{(l)}$ to $h_i^{(l+1)}$, $i$'s embedding becomes richer in a sense discussed below.

Connections between neurons in different layers are determined by $\bm{A}$ through the “message-passing” architecture

where $\Phi_{0l}(\cdot),\Phi_{1l}(\cdot)$ are parameterized functions with codomain $\mathbb{R}^H$, examples of which are provided below. We initialize $h_i^{(0)} = X_i$, interpreted as the naive node embedding that incorporates no network information. These are the usual controls under SUTVA. In subsequent layers, unit $i$'s node embedding is updated as a function of its neighbors' embeddings in the previous layer and therefore incorporates increasingly more network information as $l$ increases. This is a modeling strategy common to modern neural network architectures, that of representing a complex object in the problem domain as a Euclidean vector with learnable parameters arranged hierarchically, be it a node's network position (GNNs), a subregion of an image (CNNs), or the meaning of a word or sentence (transformers).

The final hidden layer $h_i^{(L)}$ is passed to an output layer that aggregates the vector into a scalar prediction $\Phi_o(h_i^{(L)})$ for some parameterized function $\Phi_o\colon \mathbb{R}^H \rightarrow \mathbb{R}$. The GNN ultimately returns the vector of predictions $(\Phi_o(h_i^{(L)}))_{i=1}^n \in \mathbb{R}^n$. The GNN estimator, formally defined below, amounts to a regression of $Y_i$ on $i$'s node embedding $h_i^{(L)} \in \mathbb{R}^H$, which itself depends on unknown parameters through $\{(\Phi_{0l}, \Phi_{1l})\}_{l=1}^L$.

The choices of $\Phi_{0l}(\cdot)$ and especially $\Phi_{1l}(\cdot)$ define different GNN architectures, two of which we discuss next. First we need to introduce the multilayer perceptron (MLP), which defines a sieve space commonly used in the literature. The GNN output layer $\Phi_o(\cdot)$ is typically an MLP, and $\Phi_{0l}(\cdot),\Phi_{1l}(\cdot)$ are commonly defined using MLPs.

GNN Estimator

To construct a GNN estimator, the econometrician must specify $L$, $\{(\Phi_{0l},\Phi_{1l})\}_{l=1}^L$, and $\Phi_o(\cdot)$. The latter is often taken to be an MLP ((ref)). For the hidden layers, the GIN architecture in (ref) uses sum aggregation and has a nonparametric justification, while (ref) generalizes GINs to allow for multiple aggregation functions. We use the latter in our numerical illustrations. (ref) in (ref) provides suggestions for choosing $L$.

Let $\mathcal{F}_{\text{GNN}}(L)$ denote the set of all GNNs with $L$ layers ranging over all possible functions $\{(\Phi_{0l},\Phi_{1l})\}_{l=1}^L$ and $\Phi_o(\cdot)$ in the chosen classes. Letting $d_{kl}$ be the number of parameters in $\Phi_{kl}(\cdot)$ for $k \in \{0,1\}$ and $d_o$ the number of parameters in $\Phi_o(\cdot)$, the number of parameters in a given $F \in \mathcal{F}_\text{GNN}(L)$ is $d_o + \sum_{l=1}^L (d_{0l}+d_{1l})$. Recall that each such $F$ takes as input $(\bm{X}, \bm{A})$ and outputs a vector of unit-level scalar predictions $\mathbb{R}^n$. We let $F(i,\bm{X},\bm{A}) = \Phi_o(h_i^{(L)})$ denote the $i$th component of $F(\bm{X}, \bm{A})$. A {\em GNN estimator} is a function in this set that minimizes a loss function $\ell(\cdot)$:

When $Y_i$ is real-valued, we may use squared-error loss

in which case $\hat{F}_{\text{GNN}}(\bm{X},\bm{A})$ estimates $F^*(\bm{X},\bm{A}) = ({\bf E}[Y_i \mid \bm{X}, \bm{A}])_{i=1}^n$. For binary $Y_i$, we may use the logistic loss

in which case $\hat{F}_{\text{GNN}}(\bm{X},\bm{A})$ estimates the log odds $F^*(\bm{X},\bm{A}) = (\log ({\bf E}[Y_i \mid \bm{X}, \bm{A}]/(1-{\bf E}[Y_i \mid \bm{X}, \bm{A}])))_{i=1}^n$.

Returning to the doubly robust estimator in (ref), to estimate the outcome regression with $\mathbb{R}$-valued outcomes, we restrict the sum in (ref) to the set of units $i\in\mathcal{M}_n$ for which $T_i=t$ and use squared-error loss to obtain $\hat{\mu}_t(i,\bm{X},\bm{A}) = \hat{F}_{\text{GNN}}(i,\bm{X},\bm{A})$. To estimate the generalized propensity score, we replace $Y_i$ in (ref) with $\bm{1}\{T_i=t\}$ and use logistic loss to obtain

Invariance

Modern neural network architectures often incorporate prior information in the form of input symmetries to reduce the dimensionality of the parameter space bronstein2021geometric. Convolutional neural networks (CNNs), widely used in image recognition, process grid-structured inputs and impose translation invariance. GNNs process graph-structured inputs and impose permutation invariance, as previously discussed. We next explain the purpose of this shape restriction and its economic content.

Define a {\em permutation} as an bijection $\pi\colon \mathcal{N}_n \rightarrow \mathcal{N}_n$. Abusing notation, we write $\pi(\bm{X}) = (X_{\pi(i)})_{i=1}^n$, which permutes the rows of matrix $\bm{X}$ according to $\pi$, and similarly define $\pi(\bm{D})$ and permutations of other such arrays. Likewise, we write $\pi(\bm{A}) = (A_{\pi(i)\pi(j)})_{i,j\in\mathcal{N}_n}$, which permutes the rows and columns of the matrix $\bm{A}$. We say a function $F(i,\bm{X},\bm{A})$ is {\em (permutation-)invariant} if for any bijection $\pi\colon \mathcal{N}_n\rightarrow\mathcal{N}_n$ and input $(i,\bm{X},\bm{A})$, we have $F(i,\bm{X},\bm{A}) = F(\pi(i),\pi(\bm{X}),\pi(\bm{A}))$.

Invariance appears to be necessary for estimation because the nuisance functions would otherwise depend on the unit label $i$. Computing $\hat\tau(t,t')$ would then require estimating $n$ distinct propensity scores $(p_t(i, \bm{X}, \bm{A}))_{i=1}^n$, which is infeasible when the data consists of $n$ units. The literature on neighborhood interference (ref) avoids this problem by imposing the additional restriction that $p_t(i, \bm{X}, \bm{A}) = p(W_i)$ for some function $p(\cdot)$ that does not depend on $i$. That is, units are observationally equivalent if they have identical controls $W_i$, and observationally equivalent units have identical probabilities of being assigned to an exposure mapping realization of $t$.

The analogous requirement for our setting is that assignment probabilities are equivalent for units that have isomorphic network positions, which is just another way of saying that the propensity score is invariant. This is weaker than the restriction employed under neighborhood interference.

If the propensity score is invariant, this reduces our problem from estimating $n$ separate scores to estimating only one because, for any $i$, there exists a permutation $\pi_i$ (in particular the one that only permutes units $1$ and $i$) such that $p_t(i, \bm{X}, \bm{A}) = p_t(1, \pi_i(\bm{X}), \pi_i(\bm{A}))$ and similarly for $\mu_t(\cdot)$. The right-hand side is a function $p_t(1, \cdot, \cdot)$ that does not depend on $i$, so evaluating $i$'s propensity score is now only a matter of supplying the correct $i$-specific inputs $(\pi_i(\bm{X}), \pi_i(\bm{A}))$.

The next result demonstrates that invariance is a weak requirement holding under minimal exchangeability conditions on the structural primitives.

The first three equations impose invariance on $(f_n,g_n,h_n)$, meaning that the functions do not depend directly on unit labels. This is the case for all exposure mappings $f_n$ used in the literature. Applied to $(g_n,h_n)$, invariance only says that unit identities do not influence behavior beyond the model primitives $(\bm{A}, \bm{X}, \bm{\varepsilon}, \bm{\nu})$. The final requirement says that these model primitives are themselves distributionally invariant, which is a weak condition since unit labels carry no intrinsic meaning in our setting. Distributional invariance fails to hold in design-based settings since these condition upon $\bm{\varepsilon}$, among other primitives. Then $i$'s outcome regression depends on $\varepsilon_i$, so it is not invariant with respect to observables alone.

Asymptotic Theory

We characterize the asymptotic properties of the estimators in (ref) under a sequence of models and GNNs sending $n\rightarrow\infty$. Along this sequence, the following objects may change, subject to the conditions imposed below: the exposure mapping $f_n$ from (ref), the structural functions $(g_n,h_n)$ from (ref), the distribution of the model primitives $(\bm{A}, \bm{X}, \bm{\varepsilon}, \bm{\nu})$, and the GNN functions $\{(\Phi_o, \Phi_{0l}, \Phi_{1l})\}_{l=1}^L$ and number of layers $L$ defined in (ref). Define

whose average over $i\in\mathcal{M}_n$ is the doubly robust moment condition, and let

Part (b) requires sufficient overlap for the propensity scores. Under (ref), this typically holds if $\bm{\Gamma}$ is a bounded set. Part (b) further imposes overlap on the propensity score estimator, which is common in the double machine learning literature chernozhukov2018double,farrell2015robust,farrell2021deep. It also requires that $\mathcal{M}_n$ is a nontrivial subset of $\mathcal{N}_n$. Part (c) is a standard non-degeneracy condition.

Let $Z_i = (i, Y_i, T_i, \bm{X}, \bm{A})$, and for any $x\in\mathbb{R}$, define the functions $\Psi_{p_t}(Z_i, x) = \bm{1}\{T_i=t\} (Y_i - \mu_t(i,\bm{X},\bm{A}))(x^{-1} - p_t(i,\bm{X},\bm{A})^{-1})$ and $\Psi_{\mu_t}(Z_i, x) = (x - \mu_t(i,\bm{X},\bm{A})) (1 - \bm{1}\{T_i=t\} p_t(i,\bm{X},\bm{A})^{-1})$, both of which are mean zero for any $x$.

Part (a) requires GNNs to converge at a faster than $n^{-1/4}$ rate, while part (b) is a “stochastic equicontinuity” condition. These are familiar requirements in the i.i.d.\ setting farrell2015robust. Primitive conditions are beyond the scope of the existing literature for reasons discussed in (ref), but we provide two results that may prove useful for verifying (ref). In (ref), we show that the nuisance functions can be well-approximated by $L$-neighborhood analogs for $L=O(\log n)$, making estimation a low-dimensional problem. In (ref), we reframe and combine several theoretical results in the GNN literature to show that the nuisance functions can be approximated by GNNs if they lie in a certain nonparametric subclass of invariant functions.

The next assumption is used to show that $\{\varphi_{t,t'}(i)\}_{i=1}^n$ is $\psi$-dependent kojevnikov2021limit, which enables the application of a CLT. Define

respectively $i$'s $s$-neighborhood boundary and the $k$th moment of the $s$-neighborhood boundary size. Let

where $\Lambda_n(i,s/4)$ is a constant defined in the next assumption. The second quantity is roughly a measure of the density of $s$-neighborhoods. The third quantity is a measure of dependence in that it bounds the covariance between $\varphi_{t,t'}(i)$ and $\varphi_{t,t'}(j)$ when $\ell_{\bm{A}}(i,j) \leq s$. Lastly, define

Parts (a) and (b) are used to establish that $\{\varphi_{t,t'}(i)\}_{i=1}^n$ is $\psi$-dependent. This is a nontrivial task because $(Y_i,T_i)$ are complex functions of treatments which are complex functions of unobservables. Part (b) is a Lipschitz condition that holds if potential outcomes are uniformly bounded. In particular we can take $\Lambda_n(s) = 2M$ where $M$ is the uniform bound on the ranges of $\{g_n\}_{n\in\mathbb{N}}$. Part (a) can be relaxed to allow for some forms of $\psi$-dependent unobservables under additional smoothness conditions on $g_n$ and $h_n$ kojevnikov2021limit. The simulation study in (ref) provides evidence that our methods can handle network-dependent unobservables.

The key requirement is part (d). This corresponds to Condition ND of kojevnikov2021limit, which they utilize to establish a CLT. It requires covariances to decay with network distance (controlled by $\psi_n(s)$) faster than network neighborhoods densify (controlled by $c_n(s, m; k)$). This is analogous to spatial CLTs that require mixing coefficients to decay with spatial distance $s$ faster than a particular polynomial rate, which corresponds to the rate at which $s$-neighborhoods fill with units in Euclidean space. We verify (ref) for different classes of graphs in (ref).

The next assumption is used to characterize the asymptotic properties of $\hat\sigma^2$.

Part (a) strengthens (ref)(a) to uniformly bounded outcomes. Part (b) is Assumption 4.1(iii) of kojevnikov2021limit. Part (c) is used to handle the estimation error in the nuisance functions. Under spatial dependence, one can derive an explicit upper bound on the rate at which the HAC bandwidth $b_n$ must diverge because $b_n$-neighborhoods grow at a known rate $b_n^d$ in $\mathbb{R}^d$. This is not the case for networks; (ref)(b)--(c) only implicitly characterize the rate for $b_n$ as a function of the network topology, in particular the sizes of network $b_n$-neighborhoods. It then remains to derive an explicit bandwidth that satisfies these conditions for reasonable classes of graphs. In (ref), we show that (ref) does.

Receptive Field

This section provides conditions under which the nuisance functions have low-dimensional structure that makes estimation feasible and justifies the use of shallow GNN architectures. As discussed in (ref), the number of layers $L$ in a GNN determines its receptive field $(\bm{X}_{\mathcal{N}(i,L)}, \bm{A}_{\mathcal{N}(i,L)})$ used to construct $i$'s estimate. In practice, it is common to choose a small value, which results in a receptive field that excludes most of the network. This has been found empirically to achieve better predictive performance than deep architectures with large $L$, as we will discuss in (ref). We next provide conditions under which the the GNN “estimand” $p_t(i,\bm{X}_{\mathcal{N}(i,L)}, \bm{A}_{\mathcal{N}(i,L)})$ can approximate the target $p_t(i,\bm{X},\bm{A})$ for relatively small values of $L$.

Consider the sequence of models and GNNs in (ref), and let $L=L_n$ be a possibly diverging sequence of GNN depths. Recall from (ref) that a GNN has $d_o + \sum_{l=1}^L (d_{0l}+d_{1l})$ parameters. We provide sufficient conditions for the following properties to hold for some $L = O(\log n)$.

In (ref), we discuss how these properties aid in verifying (ref). We next provide intuition by drawing an analogy to approximate sparsity conditions in the lasso literature.

Approximate Sparsity

Let $h(X_i) = {\bf E}[Y_i \mid X_i]$, and consider a lasso regression of $Y_i$ on a vector of basis functions $P(X_i)$. For the lasso prediction $P(X_i)'\hat\beta$ to be a good estimate of $h(X_i)$, we require

To verify this, it is common to impose approximate sparsity, which consists of the following two conditions belloni2014inference.

That is, $h(\cdot)$ can be approximated by a simpler function $P(\cdot)'\beta$, and estimating the latter is a low-dimensional problem. These conditions mirror (ref) and (ref) and simplify the task of showing (ref) to establishing $n^{-1} \sum_{i=1}^n (P(X_i)'(\hat\beta - \beta))^2 = o_p(n^{-1/2})$ for low-dimensional $\beta$.

The main idea in our setting is that the dependence of $Y_i$ and $D_i$ on other units decays with network distance under ANI ((ref)). That is, these quantities are primarily functions of $(\bm{X}_{\mathcal{N}(i,L)}, \bm{A}_{\mathcal{N}(i,L)})$ for some small radius $L$, which is analogous to (ref). We may then approximate $p_t(i,\bm{X},\bm{A})$ with the lower-dimensional estimand $p_t(i,\bm{X}_{\mathcal{N}(i,L)},\bm{A}_{\mathcal{N}(i,L)})$, which an $L$-layer GNN directly estimates.

Primitive Conditions

For $p_t(i,\bm{X},\bm{A})$ to approximate $p_t(i,\bm{X}_{\mathcal{N}(i,L)},\bm{A}_{\mathcal{N}(i,L)})$, we also require a form of conditional independence since the latter drops from the conditioning event the subnetwork external to the $L$-neighborhood. While ANI says that outcomes and treatments are mainly functions of units in a relatively small neighborhood, their primitives may be correlated with those of units far from the neighborhood. We consider it reasonable to assume, in the spirit of ANI, that this correlation instead decays with distance, which is the substance of the next condition.

When $r_\lambda$ is the identity function, the assumption requires that the unobservables of an $s$-neighborhood are approximately conditionally independent of the network outside of this neighborhood, where the approximation error is shrinking with the radius $s$. More generally, it allows the $s$-neighborhood to be approximately conditionally independent of the network outside the greater $r_\lambda(s)$-neighborhood for $r_\lambda(s) \geq s$.

The result allows $L$ to be of order $\log n$, which quantifies the sense in which the GNN can be shallow. The lower bounds on $L$ given in the theorem are not feasible choices, being dependent on unknowns $r_\lambda(\cdot)$ and $\alpha$. This is similar to how finite-sample bounds for the lasso require restrictions on the penalty parameter involving unknown constants. In (ref), we illustrate the performance of different choices of $L$ in simulations and provide suggestions for choosing $L$ in practice.

The first half of the theorem concerns the propensity score. Exponential decay of the interference bounds in (ref) holds for the models discussed in (ref). Real-world networks are typically sparse, usually formalized as $n^{-1} \sum_{i=1}^n n(i,1) = O_p(1)$, which the theorem mildly strengthens to a second-moment condition. The second half of the theorem concerns $\mu_t(\cdot)$. The requirement $n^{-1} \sum_{i=1}^n n(i,s)^2 = O_p(e^{\xi s})$ for $\xi<2\alpha$ says $s$-neighborhoods grow at a rate $\xi$ not too much larger than the rate $\alpha$ at which interference decays. A similar type of condition is required for a central limit theorem, as discussed in (ref) and (ref).

The last result provides sufficient conditions for (ref).

Simulation Study

We design a simulation study with three objectives in mind. The first is to illustrate the finite-sample properties of our proposed estimators for different choices of $L$. The second is to demonstrate that shallow GNNs can perform well even on “wide” networks that otherwise may require many layers in the absence of the low-dimensional structure established in (ref). The third is to compare the performance of GNNs to that of standard machine learners utilizing prespecified controls based on (ref).

Design

We simulate $\bm{A}$ from two random graph models. The random geometric graph model sets $A_{ij} = \bm{1}\{\lVert\rho_i-\rho_j\rVert \leq r_n\}$ for $\{\rho_i\}_{i=1}^n \stackrel{iid} \sim \mathcal{U}([0,1]^2)$ and $r_n = (5/(\pi n))^{1/2}$, where $\pi$ is the transcendental number. The Erd\H{o}s-R\'{e}nyi model sets $A_{ij} \stackrel{iid}\sim \text{Ber}(5 / n)$. Both have limiting average degree equal to five. The former model results in “wide” networks with high average path lengths that grow at a polynomial rate with $n$, while the latter results in low average path lengths of $\log n$ order. For $n=2000$, the average path length is about 39.5 for random geometric graphs and 4.9 for Erd\H{o}s-R\'{e}nyi graphs.

We draw mutually independent sequences $\{X_i\}_{i=1}^n \stackrel{iid}\sim \mathcal{U}(\{0,0.25,0.5,0.75,1\})$, $\{\varepsilon_i\}_{i=1}^n \stackrel{iid}\sim \mathcal{N}(0,1)$, and $\{\nu_i\}_{i=1}^n \stackrel{iid}\sim \mathcal{N}(0,1)$ independently of $\bm{A}$. For some vectors $\bm{D}=(D_i)_{i=1}^n$ and $\bm{\nu}=(\nu_i)_{i=1}^n$, define

where $\theta = (\alpha,\beta,\delta,\gamma)$. We generate $\{Y_i\}_{i=1}^n$ from a linear-in-means model with $Y_i=V_i(\bm{Y},\bm{\varepsilon};\theta_y)$ and $\theta_y=(0.5,0.8,10,-1)$. We generate $\{D_i\}_{i=1}^n$ according to (ref), so that $D_i = \bm{1}\{V_i(\bm{D},\bm{\nu};\theta_d) > 0\}$ with $\theta_d=(-0.5,1.5,1,-1)$. The equilibrium selection mechanism is myopic best-response dynamics starting from the initial condition $\{D_i^0\}_{i=1}^n$ for $D_i^0 = \bm{1}\{V_i(\bm{0},\bm{\nu};\theta_d) > 0\}$.

The design induces a greater degree of dependence than what our assumptions allow. The error term $\nu_i + \sum_{j=1}^n A_{ij} \nu_j / \sum_{j=1}^n A_{ij}$ is not conditionally independent across units unlike what (ref)(a) requires. Also, back-of-the-envelope calculations indicate that peer effects are sufficiently large in magnitude that (ref)(d) is violated. Section (ref) presents results for a design that satisfies our assumptions.

The estimand is $\tau(1,0)$ in (ref), whose true value is zero. About 57 percent of units select into treatment, so the effective sample size used to estimate the outcome regressions is around $n/2$ since $\mu_t(i,\bm{X}, \bm{A})$ is estimated only with observations for which $T_i=t$. We report results for $n=1000,2000,4000$.

Nuisance Functions

For the GNNs, we use the PNA architecture in (ref) with aggregator $\Gamma_2(\cdot)$ defined in the example and $L=1,2,3$. Both $\phi_{0l}$ and $\phi_{1l}$ are one-layer MLPs ((ref)) with widths $H=3,6,9$. We optimize the GNNs with full-batch gradient descent using the default PyTorch implementation of the “Adam” gradient descent algorithm NEURIPS2019_9015. We randomize initial parameter values and utilize a learning rate of 0.01. For $\phi_{1l}(\cdot)$, we use a linear layer ((ref)), which is the default for the {\tt PNAConv} class in the PyTorch Geometric package fey2019fast. For $\phi_{0l}(\cdot)$, we use a one-layer MLP with ReLU activation. Finally, $\Phi_o(\cdot)$ is a linear layer.

We compare GNNs to conventional machine learners using the prespecified controls given in (ref), which are analogous to those used in the numerical illustrations of emmenegger2022treatment and forastiere2021identification. We report results for two different machine learners. The first is an MLP with the same number of layers $L$ and width $H$ as the GNNs and trained using the same gradient descent algorithm and learning rate. The second is a random forest with 500 trees and minimum leaf size 10.

Results

Tables (ref) and (ref) report the results of 5000 simulations for the random geometric graph and Erd\H{o}s-R\'{e}nyi models, respectively. Rows marked $\hat\tau$, SE, and CI report the average point estimate, average standard error, and empirical coverage of 95-percent confidence intervals, respectively. Since the true estimand is zero, the absolute value of $\hat\tau$ also equals the bias. Rows SE+ and CI+ use the gao2025causal standard errors.\footnote{In the notation of (ref), their variance estimator is $m_n^{-1}\tilde{\bm{\tau}}'\bm{K}_+^U\tilde{\bm{\tau}}$ where $\bm{K}_+^U = \bm{Q} \max\{\bm{\Lambda},\bm{0}\} \bm{Q}'$ and $\bm{Q}\bm{\Lambda}\bm{Q}'$ is the eigendecomposition of $\bm{K}^U$.} The “oracle” rows report results using the true standard error, computed by taking the standard deviation of $\hat\tau(1,0)$ across simulation draws. The “naive” rows report results using i.i.d.\ standard errors, which illustrate the degree of dependence. All other SE and CI rows use our variance estimator in (ref).

We first compare the standard machine learners to GNNs with $L=1$ since all share the same receptive field in this case. The bias obtained with GNNs is smaller, often by half, particularly in the Erd\H{o}s-R\'{e}nyi case. The most competitive alternative is random forests which achieves a similar bias when the network is the random geometric graph, but in the case of the Erd\H{o}s-R\'{e}nyi model, its bias is over 30 percent larger. This suggests that GNNs can learn a different function of $(\bm{X},\bm{A})$ than common prespecified controls, one that can better adjust for confounding. The improvement in bias using GNNs does not come at an apparent cost to variance as the SEs are comparable across methods.

Second, we compare the GNN estimators across different choices of $L$. The best performance in terms of bias is achieved with $L=2$. This is the case for both random graph models and is particularly notable for the random geometric graph because its average path length is substantially larger than the radius of the receptive field. This demonstrates that GNNs can perform well despite only controlling for $(\bm{X}_{\mathcal{N}(i,2)}, \bm{A}_{\mathcal{N}(i,2)})$, which is possible due to the low-dimensional structure established in (ref). Unsurprisingly, $L=2$ outperforms $L=1$ since the latter only adjusts for 1-neighborhood confounding. In principle, $L=3$ accounts for higher-order network confounds, but the bias turns out to be slightly larger and the coverage slightly worse.

The choice of $L=2$ is not unusual in the GNN literature. For example, das2023credit and ying2018graph use $L=2$, while ma2021graph use $L=3$. zhou2021understanding compute the prediction error of GNNs across several different datasets and architectures with $L=2,4,8,\dots$ and find that $L=2$ has the best performance. The fact that GNN performance often fails to improve, and indeed can worsen, with larger $L$ is well known in the GNN literature.\footnote{bronstein2020do writes, “Significant efforts have recently been dedicated to coping with the problem of depth in graph neural networks, in hope to [sic] achieve better performance and perhaps avoid embarrassment in using the term `deep learning' when referring to graph neural networks with just two layers.”} See (ref) for a survey of possible explanations.

The oracle CIs achieve coverage close to the nominal level across most sample sizes and architectures, which illustrates the quality of the normal approximation. Our SEs result in some undercoverage, which is common for HAC estimators, though coverage tends toward the nominal level as $n$ grows for $L=2$. The gao2025causal SEs are generally larger than the oracle SEs and hence conservative, which improves coverage in smaller samples.

Empirical Application

We revisit the analysis of he2024measuring on the diffusion of microfinance through rural villages in Karnataka, India. They use data collected by banerjee2013diffusion containing twelve dimensions of social relationships, demographic details, and microfinance adoption decisions from 43 villages involved with Bharatha Swamukti Samsthe's (BSS) microfinance program in 2007. BSS initiated the program by meeting with a select group of village “leaders” who were asked to spread the word about microfinance.

Following the analysis of he2024measuring, the unit of observation is the household, and household observables $X_i$ include the normalized total number of households within each village and indicators for participation in self-help groups, savings activities, and caste composition. They construct three social networks from the multigraph data: $G_{ee}$ represents connections through material exchanges like borrowing or lending essentials; $G_{sc}$ captures social activities including advice sharing or joint religious attendance; and $G_{all}$ is the union of $G_{ee}$ and $G_{sc}$. We report results with $\bm{A}$ set to each of the three. The average path length is about 3 for all graphs.

We consider three different definitions of the treatment. In the “leader case,” $D_i$ is an indicator for whether household $i$ has a leader. In the “leader-adopter case,” it is an indicator for whether the household has a leader who adopts microfinance in the first trimester of the study. These definitions are used by he2024measuring. We add to these the “adopter case,” where $D_i$ is simply an indicator for whether any household member adopts microfinance in the first trimester. The outcome $Y_i$ is an indicator for whether an individual in household $i$ adopts microfinance starting in the first trimester of the study or later.\footnote{In the adopter case, the treatment time period intersects with that of the outcome in the first trimester, which generates unwanted feedback between the outcome and selection models. To avoid this, we should define $Y_i$ as an indicator for adoption {\em after} the first trimester. Fortunately, adoption decisions are never reversed in the data, so this coincides with the present definition of $Y_i$.}

Comparison with He and Song (2024)

We next discuss estimands in the context of the adopter case. he2024measuring propose a novel estimand called “average diffusion at the margin” (ADM). This measures the expected number of neighbors induced to adopt microfinance as a result of the ego's adoption. To identify the ADM, they assume the following. First, initial adoption decisions are unconfounded (their Assumption 2.1), and in the application, they only use household observables $X_i$ as the controls. Second, their selection model, as described in (ref), is a parametric single-agent discrete choice model. Third, they assume adoption decisions are irreversible in that $Y_i \geq D_i$, which is true in the application. Finally, as discussed in (ref), the econometrician must specify the maximum number of rounds of diffusion that take place between the measurement of $D_i$ and $Y_i$, and they choose $K=1$ in the application.

Our approach has several advantages. First, we do not require knowledge of the number of within-period rounds of diffusion, binary outcomes, or irreversibility of adoption decisions. Second, we use a richer set of network controls that includes covariates of higher-order neighbors, without assuming a known function $W_i$ as in (ref). Third, we employ a nonparametric selection model allowing for peer effects in initial adoption. Fourth, our outcome model (ref) also allows for peer effects in subsequent adoption, as well as higher-order diffusion beyond the $K=1$ neighborhood since outcomes depend on the entire initial adoption vector $\bm{D}$. ANI posits that this dependence decays with distance, which is a feature of most diffusion models, as information from distant units is less likely to diffuse to the ego.

The cost of imposing less structure than he2024measuring is that the ADM may not be identified under our assumptions. Instead, we consider two estimands defined by the following exposure mapping:

with $\mathcal{M}_n$ defined in (ref). Then $\tau(1,0)$ ($\tau(2,0)$) measures the effect on own adoption of going from 0 to 1 (two or more) adopting neighbor(s) for units with at least one neighbor.\footnote{Assumption 2.1 of he2024measuring implies that $\{D_i\}_{i=1}^n$ is independent conditional on observables. Our estimands have a causal interpretation under the same condition by (ref).} This sheds light on a different dimension of diffusion relative to the ADM. Whereas the ADM measures how many alters are affected by the ego's adoption, our estimands quantify the effect of having multiple adopting neighbors on the ego's adoption. We will find that having multiple adopting neighbors has a much larger effect than having only one.

As previously stated, he2024measuring define the treatment in two ways. One is a binary indicator for having a leader in the household, the idea being that all leaders were initially informed about microfinance and told to spread the word. However, not all leaders adopted in the first trimester, which perhaps motivates the second definition, a binary indicator for having an adopting leader in the household. In our view, it may be plausible to argue that microfinance adoption in the initial period is as good as random within observable subpopulations, but it is less plausible to make the same case for being a leader, which is likely determined by a complex social process. We therefore consider a third definition, which is simply an indicator for adopting microfinance in the initial period, irrespective of having a leader in the household. Recall that the interpretations of the causal estimands above pertain to the third definition; in our view, the interpretations are less clear when treatment is defined as in the other cases.

Results

We present estimates of $\tau(1,0)$, $\tau(2,0)$, and the ADM for the three network specifications and three treatment definitions introduced above. We use three different machine learning estimators of the nuisance functions. The first is GNNs using the same PNA architecture, learning rate, and gradient descent algorithm as the simulation study (see (ref)). The other estimators use the same prespecified controls (ref) and machine learners (MLPs and random forests) as the simulation study. For both neural networks, we set the width to $H=4$ to match the number of household covariates $X_i$.

To compute the estimates, we concatenate the village networks into a single adjacency matrix of size $n=4413$. For the non-ADM estimates, we trim observations with propensity scores outside of $[0.01,0.99]$.\footnote{The only substantial trimming occurs in the leader-adopter case for $\tau(2,0)$. This is because the number of units with two or more treated units is small, no more than 100. The worst case is network $G_{ee}$ for which the smallest post-trimming sample size across $L$ is 1020 for GNNs, 889 for MLPs and 554 for random forests. Fortunately, these are sufficiently large that we draw similar conclusions from the estimates regardless of $L$ or the network.} Standard errors are obtained from the variance estimator defined in (ref). Across all network definitions, the bandwidth (ref) equals one.

Tables (ref) and (ref) report results for $\tau(1,0)$ and $\tau(2,0)$ respectively, as well as ADM estimates. First consider $\tau(1,0)$, which contrasts microfinance adoption rates for units with 1 versus 0 initially adopting neighbors. The GNN results are similar across $L$. For the leader case, we obtain precise zeros for almost all estimates. For the leader-adopter case, the GNN estimates are substantially smaller in magnitude than those of the ADM, at most half the magnitude. While the ADM is a different estimand, $\tau(1,0)$ is perhaps the case where they are most logically comparable, and the smaller effect sizes we find may be due to the use of richer network controls. For the adopter case, the GNN estimates are an order of magnitude smaller than the ADM estimates. The MLP estimates are comparable to the GNN estimates, while the random forest estimates are generally smaller, particularly in the adopter case.

The estimand using $\tau(2,0)$ contrasts units with $2+$ versus 0 initially adopting neighbors. The estimates for the leader case are precise zeros, but unlike $\tau(1,0)$, we find sizeable effects for the leader-adopter and adopter cases. For the latter, the estimates are around 20 percentage points. The MLP estimates are comparable to the GNN estimates, while the random forest estimates are substantially smaller. Regardless of method, the conclusion from both tables is that the effect of having multiple adopting neighbors is more than triple the effect of having only one, but the magnitudes are less than half the ADM estimates.

Conclusion

Existing work on network interference with unconfoundedness assumes that it suffices to control for a known, low-dimensional function of the network and covariates $W_i$. In this respect, the approaches may be viewed as semiparametric. We propose to use GNNs to nonparametrically learn the function and provide a behavioral model under which it is low-dimensional and estimable with shallow GNNs.

Our contributions are twofold. First, we observe that the standard formulation of unconfoundedness that controls for a known function of the network is limited in terms of microfoundations, ruling out simultaneous peer effects in the outcome and treatment selection stages. We consider a model of approximate neighborhood interference (ANI) that allows for both. The challenging aspect of the model is that it induces high-dimensional network confounding, which is presumably why it has not been previously studied in a fully nonparametric setting.

Second, we make the case that estimation remains feasible. We draw an analogy between ANI and approximate sparsity conditions in the lasso literature, which posit that a high-dimensional regression function is well-approximated by a lower-dimensional analog. Under ANI and additional conditions, we show that the propensity score and outcome regression can be approximated by low-dimensional functions of the ego's $L$-neighborhood network for relatively small $L$. leung2022causal studies the implications of ANI for asymptotic inference in randomized control trials, while we highlight its utility for handling high-dimensional network confounding.

\part{Supplementary Appendix}

\makeatletter \@addtoreset{section}{part} \makeatother \setcounter{section}{0} \numberwithin{equation}{section}

Theoretical Properties of GNNs

(ref)

{\bf Mean-Squared Error.} This section discusses verification of (ref), which is currently beyond the scope of the literature. In the notation of the (ref), let $\varrho$ denote either the propensity score $p_t$ or outcome regression $\mu_t$. To verify (ref)(a), (ref) implies that it suffices to show

for $L = O(\log n)$. This should be more feasible to verify directly given that (a) $\hat{\varrho}$ is an $L$-layer GNN which only uses information from $(\bm{X}_{\mathcal{N}(i,L)}, \bm{A}_{\mathcal{N}(i,L)})$, (b) the GNN is low-dimensional under the conditions of (ref), (c) $\varrho$ is invariant like $\varrho$ under the conditions of (ref), and (d) the data $\{(Y_i, T_i, \bm{X}_{\mathcal{N}(i,L)}, \bm{A}_{\mathcal{N}(i,L)})\}_{i=1}^n$ is $\psi$-dependent conditional on $(\bm{X}, \bm{A})$ by an argument similar to (ref).

Recent work by wang2024graph provides sufficient conditions for (ref) (see their section 3.2), but some limitations prevent a direct application to our setup.

In (ref) below, we provide a characterization, drawing heavily from existing theory in the GNN literature. The result says that, for any fixed $n$ and $L$, if the function lies within a certain nonparametric subclass of invariant functions, then there exists a sequence of GNNs converging to it. This is a step toward what is eventually needed, which is a rate of convergence in terms of the GNN parameters wang2024graph.

{\bf Stochastic Equicontinuity.} (ref)(b) is typically established from (ref)(a) and additional conditions. It is most straightforward to verify with cross-fitting, but it can be established without cross-fitting for certain machine learners in low-dimensional settings. farrell2021deep do so for MLPs and i.i.d.\ data (their Lemma 10). wang2024graph provide sufficient conditions for GNNs (their Appendix D).

In our formulation, we explicitly state the stochastic equicontinuity (SE) condition for the propensity score. This is often left implicit in standard i.i.d.\ settings, but it requires careful attention here. Specifically, the condition for $\Psi_{\mu_t}$ is standard and concerns the outcome regression, while the condition for $\Psi_{p_t}$ is the analog for the propensity score. The latter is not explicitly stated in the literature because, when the data is i.i.d.\ and SUTVA holds, it is straightforward to verify from first principles farrell2015robust. However, due to the complexity of our setting, verification from first principles is not apparently possible without stronger restrictions.

To see this, we discuss two methods of verifying SE for $\Psi_{p_t}$. The first is to observe that $\Psi_{\mu_t}$ and $\Psi_{p_t}$ have the same structure, so techniques used by wang2024graph to verify SE for the former can be applied to the latter. The second is to either use cross-fitting or restrict interference in the outcome model to enable verification from first principles.

{\bf Method 1.} Appendix D of wang2024graph verifies SE for $\Psi_{\mu_t}$. Under their assumptions, the arguments are directly applicable to $\Psi_{p_t}$ because they only rely on the following properties that are shared by both functions given uniformly bounded outcomes.

{\bf Method 2.} The next two lemmas show that the SE condition for $\Psi_{p_t}$ holds if either the GNNs are trained using cross-fitting or outcomes follow a generalized neighborhood interference model (ruling out endogenous peer effects in the outcome but not the selection stage).\footnote{By comparison, Corollary 4 of wang2024graph verifies the condition using independent clusters (their \S 4.2.2) and SUTVA (their \S 4.1), which is a special case of neighborhood interference. Lemma E.4 of emmenegger2022treatment verifies SE using cross-fitting (their Algorithm 1) and neighborhood interference (their equation (1)).} With a single network, it is not generally possible to cross-fit in a manner that satisfies the conditions of the next lemma, but the purpose of the exercise is to illustrate a technical point discussed in (ref) below. Without either restriction, a key step breaks down, and the challenge arises not from dependence per se but from its combination with our richer model of interference.

WL Function Class

MLPs can approximate any measurable function hornik1989multilayer, so given the discussion in (ref), a natural question is whether GNNs can approximate any measurable, {\em invariant} function of graph-structured inputs. In other words, for GNNs to approximate the propensity score or outcome regression, is it enough to assume that these functions are invariant (and satisfy appropriate regularity conditions)? For reasons related to the graph isomorphism problem, stronger restrictions appear to be necessary.

To see why, let $\mathcal{F}$ be a set of functions mapping $(\bm{X},\bm{A})$ to $\mathbb{R}$. What properties must $\mathcal{F}$ have for any invariant function to be well approximated by a sequence of functions in $\mathcal{F}$? chen2019equivalence show that there must exist $F \in \mathcal{F}$ such that $F$ can separate any pair of non-isomorphic graphs in that $F(\bm{X},\bm{A}) \neq F(\bm{X}',\bm{A}')$ for any non-isomorphic $(\bm{X},\bm{A}), (\bm{X}',\bm{A}')$. A function with this property solves the graph isomorphism problem, a problem for which no known polynomial-time solution exists kobler2012graph,morris2021weisfeiler. Since GNNs are an example of $\mathcal{F}$ that can be computed in polynomial time, this strongly suggests that approximating any invariant function is too demanding of a requirement.

To define the subclass of invariant functions that GNNs can approximate, we need to take a detour and discuss graph isomorphism tests. The subclass will be defined by a weaker graph separation criterion than solving the graph isomorphism problem, in particular one defined by the {\em Weisfeiler-Leman (WL) test}. This is a (generally imperfect) test for graph isomorphism on which almost all practical graph isomorphism solvers are based morris2021weisfeiler.

Given a labeled graph $(\bm{X},\bm{A})$, the WL test outputs a graph coloring (a vector of labels for each unit) according to the following recursive procedure, whose definition follows maron2019provably. At each iteration $t>0$, each unit $i$ is assigned a color $C_t(i)$ from some set $\Sigma$ (e.g.\ the natural numbers) according to

where $\Phi(\cdot)$ is a bijective function that takes as input a color and a multiset of neighbors' colors.\footnote{Strictly speaking, this is the 1-WL test.} Intuitively, at each iteration, two units are assigned different colors if they differ in the number of identically colored neighbors, so that at iteration $t$, colors capture some information about a unit's $(t-1)$-neighborhood. Colors are initialized at $t=0$ using a deterministic rule that assigns each $i$ to the same color $C_0(i) \in \Sigma$ if and only if they have the same covariates $X_i$. At each iteration, the number of assigned colors increases, and the algorithm converges when the coloring is the same in two adjacent iterations. This takes at most $n-1$ iterations since there cannot be more than $n$ distinctly assigned colors.

To test whether two labeled graphs are isomorphic, the procedure is run in parallel on both graphs for some number of iterations, typically until convergence. At this point, if there exists a color such that the number of units assigned that color differs in the two graphs, then the graphs are considered non-isomorphic. This procedure correctly identifies isomorphic graphs, but it is underpowered since there exist non-isomorphic graphs considered isomorphic by the WL test morris2021weisfeiler. Also, because the number of colors increases each iteration, the test is more powerful when run longer.

morris2019weisfeiler and xu2018powerful note the similarity between the GNN architecture (ref) and WL test (ref). The former may be viewed as a continuous approximation of the latter, replacing the hash function $\Phi(\cdot)$ with a learnable aggregator $\Phi_{1l}(\cdot)$. They formally show that any GNN has at most the graph separation power of the WL test and furthermore there exist architectures as powerful.

Returning to the original problem, we now define the class of functions approximated by GNNs in terms of the WL test. Let $\mathcal{S}$ denote the support of $(\bm{X},\bm{A})$.

This is essentially Definition 2 of azizian2021expressive. Intuitively, if $\mathcal{E}$ is at most as separating as $\mathcal{F}$, the latter is more complex in the sense that some function in $\mathcal{F}$ can separate weakly more elements of $\mathcal{S}$ than any function in $\mathcal{E}$.

Let $f_{\text{WL},L}$ denote the function of $(\bm{X},\bm{A})$ with codomain $\Sigma^n$ that outputs the vector of node colorings from the WL test run for $L$ iterations. Let $\mathcal{C}(\mathcal{S})$ be the set of continuous functions with domain $\mathcal{S}$. For any $L\in\mathbb{N}$, define the {\em WL function class}

This is the set of continuous functions of $(\bm{X},\bm{A})$ that are at most as separating as the WL test with $L$ iterations.

The next result says that a function can be approximated by $L$-layer GNNs under the shape restriction that they are elements of the WL function class. This is a stronger shape restriction than invariance because, by construction, the output of the WL test is invariant, so $\mathcal{F}_\text{WL}(L)$ is a subset of the set of all invariant functions.

Consider the GNN architecture in (ref) with $\phi_{0l}(\cdot), \phi_{1l}(\cdot), \Phi_o(\cdot)$ being MLPs. For technical reasons, we augment the output layer of the architecture with an additional MLP layer $L+1$ at the output stage, so for an MLP $\Phi^*\colon \mathbb{R}^n \rightarrow \mathbb{R}^n$, the GNN output is $\Phi^*(\Phi_o(h_1^{(L)}), \ldots, \Phi_o(h_n^{(L)}))$. Let $\mathcal{F}_{\text{GNN}*}(L)$ denote the set of such GNNs with $L$ layers, ranging over the parameter space of the MLPs, including their widths and depths. Finally, for any $f \in \mathcal{F}_{\text{GNN}*}(L)$, let $F(i, \bm{X}, \bm{A})$ denote the $i$th component of $F(\bm{X},\bm{A})$.

In other words, any function in the class $\mathcal{F}_\text{WL}(L)$ can be approximated by a sequence of $L$-layer GNNs in $\mathcal{F}_{\text{GNN}*}(L)$. The result is a consequence of a Stone-Weierstrauss theorem due to azizian2021expressive and a version of the morris2019weisfeiler and xu2018powerful result on the equivalent separation power of GNNs and the WL test. The proof is given below.

The result is essentially Theorem 4 of azizian2021expressive but with the distinction that they use $\mathbin{\scalebox{1.5}{\ensuremath{\cup}}}_L \mathcal{F}_{\text{GNN}*}(L)$ in place of $\mathcal{F}_{\text{GNN}*}(L)$ and $\{F_{\text{WL},\infty}\}$ in place of $\{F_{\text{WL},L}\}$. That is, their theorem states that the set of GNNs ranging over all possible numbers of layers can approximate any continuous function at most as separating as the WL test run until convergence.

(ref) states their result for fixed $L$, and the proof is straightforward from prior results. However, our framing clarifies one of the roles of depth, namely that it determines the strength of the shape restriction implicitly imposed on the function being approximated by GNNs. In particular, because the WL test is more powerful when $L$ is larger, meaning when run for more iterations, (ref) implies that deeper GNNs can approximate weakly richer function classes, or equivalently, impose weaker shape restrictions. We discuss the economic significance of this point in the next subsection.

Disadvantages of Depth

As discussed in (ref), the receptive field is the main consideration when selecting $L$, but (ref) provides a second consideration, which is imposing a weaker implicit shape restriction. It shows that, for GNNs to approximate a target function well, the target must satisfy a shape restriction stronger than invariance, namely that it is at most as separating as the WL test with $L$ iterations. The larger the choice of $L$, the weaker the shape restriction imposed. This cuts against the standard practice of choosing small values of $L$. However, there are several reasons why shallow architectures remain preferable despite this result.

{\bf Low returns to depth.} For a given graph, how many iterations are required for the WL test to converge? This corresponds to the choice of $L$ for which the shape restriction is weakest. If the number is large for most graphs in practice, then it would suggest gains to choosing large values of $L$.

Unfortunately, the answer is not generally known, being determined by the topology of the input graph in a complex manner. However, there is a range of results bounding the number of iterations required for convergence. For instance, kiefer2020iteration construct graphs for which the WL test requires $n-1$ iterations to converge, so such graphs require $n-1$ layers to obtain the weakest shape restriction. This makes the estimation problem extremely high-dimensional, requiring substantially more layers than what is typically required for the receptive field to encompass the entirety of the network.

Fortunately, theoretical and empirical evidence suggest that such examples are more the exception than the rule and that small choices of $L$ are typically enough to separate many graphs. babai1980random show that, with probability approaching one as $n\rightarrow\infty$, in an $n$-unit network drawn uniformly at random from the set of all possible networks, the WL test assigns all units different colors (recall the test must converge at this point) after only {\em two} iterations morris2021weisfeiler. Thus, roughly speaking, for large networks, the weakest possible shape restriction is generically achieved with only $L=2$. This might suggest that using a small number of layers is not restrictive in practice. Indeed, zopf20221 provide empirical evidence on this point, showing that the vast majority of graphs in their dataset can be separated using the WL test after a single iteration.

{\bf Cost of depth.} Empirically, larger $L$ has been found to result in worse predictive performance, and several explanations have been proposed. The “oversmoothing” phenomenon li2018deeper,oono2020graph posits that node embeddings tend to become indistinguishable across many units as the number of layers grows. In random geometric graphs (see (ref)), $L$-neighborhood sizes grow polynomially with $L$, while in Erd\H{o}s-R\'{e}nyi graphs, the growth rate is exponential. Accordingly, a small increase in $L$ can induce a large increase in the number of elements aggregated by $\Phi_{1l}(\cdot)$, so by a law of large numbers intuition, the resulting node embeddings tend to concentrate on the same value. Since node embeddings are meant to represent network positions, which tend to be quite heterogeneous across units, this results in poor predictive performance.

The “oversquashing” phenomenon alon2021bottleneck,topping2022understanding posits that, as $L$ grows, the GNN aggregates an exceedingly large amount of information due to the growth in neighborhood sizes. This information is stored in node embeddings of relatively small dimension $H$, resulting in information loss, so the effective size of the receptive field remains small as $L$ grows.

zhou2021understanding provide a third explanation, that certain features of common architectures are responsible for variance inflation. In fact, even weaker shape restrictions than those imposed by (ref) are possible using more complex “$k$-GNN” architectures, which would theoretically improve bias, but these have greater computational cost and empirically exhibit worse predictive performance and higher variance than the standard architecture (ref) dwivedi2022benchmarking. These disadvantages may explain in part the common use of the standard architecture with few layers.

Verifying \S8 Assumptions

leung2022causal, \S A, verifies analogs of Assumptions (ref)(d) and (ref) from an older working paper version of kojevnikov2021limit. This section repeats the exercise for the published version of the assumptions. We assume throughout that $\max\{\gamma_n(s/2),\psi_n(s)\} \leq \text{exp}(-c(1-4/p)^{-1} s)$ for some $c>0$ and $p$ in (ref)(a). As in leung2022causal, we say a sequence of networks exhibits polynomial neighborhood growth if

for some $C>0$, $d\geq 1$. The sequence exhibits exponential neighborhood growth if

for some $C > 0$ and $\beta = \log \delta(\bm{A})$ leung2022causal.

(ref)(d) and ((ref))

For polynomial neighborhood growth, choose $v_n = n^{1/(\alpha d)}$ for $\alpha > 2$. The second term in (ref) is at most $n^{3/2} \text{exp}(-c\,n^{1/(\alpha d)}) = o(1)$. The first term is at most $n^{-1/2} \sum_{s=0}^\infty (Cn^{1/\alpha}) (Cs^d) \text{exp}(-c\,s) = o(1)$ for $k=1$, and for $k=2$, it is at most $n^{-1} \sum_{s=0}^\infty (Cn^{1/\alpha})^2 (Cs^d) \text{exp}(-c\,s) = o(1)$. Finally $\eqref{3ND} \leq \sum_{s=0}^\infty Cs^d \text{exp}(-c\, s) < \infty$.

For exponential neighborhood growth, choose $v_n = \alpha\beta^{-1} \log n$, $\alpha \in (1.5\beta c^{-1}, 0.5)$, with $c$ from the definition of $\psi_n(s)$ above. Such an $\alpha$ exists only if $c>3\beta$, which requires $\psi_n(s)$ to decay sufficiently fast relative to neighborhood growth. The second term in (ref) is then at most $n^{3/2} \text{exp}(-c\alpha\beta^{-1} \log n) = n^{1.5-c \alpha\beta^{-1}} = o(1)$. For $k=1$, the first term is at most $n^{-1/2} \sum_{s=0}^\infty C^2\text{exp}(\alpha \log n) \text{exp}((\beta-c)s) = o(1)$, and for $k=2$, it is at most $n^{-1} \sum_{s=0}^\infty C^2\text{exp}(2\alpha \log n) \text{exp}((\beta-c)s) = o(1)$. Finally, $\eqref{3ND} \leq \sum_{s=0}^\infty C\text{exp}((\beta -c)s) < \infty$.

Bandwidth

We employ a mix of formal and heuristic arguments to show that the bandwidth (ref) satisfies (ref)(b)--(c). Under polynomial neighborhood growth, as argued in \S A.2 of leung2022causal, $\mathcal{L}(\bm{A}) \approx n^{1/d}$, in which case $b_n = \mathcal{L}(\bm{A})^{1/4} \approx n^{1/(4d)}$. Then (ref)(c) holds because $n^{-1} \sum_{i=1}^n n(i,b_n) = Cb_n^d \approx n^{1/4} = o(\sqrt{n})$. (ref)(b) holds because, taking $\epsilon=1-4/p$,

Under exponential neighborhood growth, as argued in \S A.2 of leung2022causal, $\mathcal{L}(\bm{A}) \approx \log n / \log \delta(\bm{A})$, in which case $b_n \approx 0.25 \log n / \log \delta(\bm{A})$. Then (ref)(c) holds because $n^{-1} \sum_{i=1}^n n(i,b_n) = C\text{exp}(\beta b_n) \approx n^{1/4}$. (ref)(b) holds because, taking $\epsilon=1-4/p$,

which is $o(1)$ if $c>\beta$, which is weaker than the requirement $c>3\beta$ in (ref).

Additional Simulation Results

(ref) presents simulation results from a modified version of the design in (ref) that exhibits sufficient weak dependence required by our assumptions. The results are very similar. In the definition of $V_i(\bm{D},\bm{\nu}; \theta)$ we replace $\nu_i + \sum_{j=1}^n A_{ij}\nu_j / \sum_{j=1}^n A_{ij}$ with a conditionally independent error term with the same distribution, namely $\nu_i + (\sum_{j=1}^n A_{ij})^{-1/2} \tilde\nu_i$ where $\{\tilde\nu_i\}_{i=1}^n \stackrel{iid}\sim \mathcal{N}(0,1)$ is independent of all other structural primitives. The graph exhibits polynomial neighborhood growth in the sense of (ref) leung2022causal. Since the data-generating process satisfies (ref) with exponential decay, Assumptions (ref) and (ref) are satisfied by the calculations in (ref).

Supporting Lemmas

Throughout this section we abbreviate $\bm{1}_i(t) = \bm{1}\{T_i=t\}$.

The following notion of weak network dependence is due to kojevnikov2021limit. For any $H,H' \subseteq \mathcal{N}_n$, define $\ell_{\bm{A}}(H,H') = \min\{\ell_{\bm{A}}(i,j)\colon i \in H, j \in H'\}$. Let $\{Z_i\}_{i=1}^n \subseteq \mathbb{R}$ be a triangular array, $\bm{Z}_H = (Z_i)_{i \in H}$, $\mathcal{L}_d$ be the set of bounded $\mathbb{R}$-valued Lipschitz functions on $\mathbb{R}^d$, $\text{Lip}(f)$ be the Lipschitz constant of $f \in \mathcal{L}_d$, and

Proofs of Main Results

\FloatBarrier \phantomsection \addcontentsline{toc}{section}{References}