Graph Neural Networks: Theory for Estimation with Application on Network Heterogeneity

Graph Neural Networks: Theory for Estimation with Application on Network Heterogeneity

Introduction

The world is inherently interconnected, whether it is the natural, engineering, or social domain. For instance, fundamental particles bond together to form larger molecules, which together form building blocks of organic or inorganic objects. Human-engineered systems, such as computers or the Internet, are made up of large amounts of small components, connected together in well-designed manners in order to perform more complex functionalities. Social animals, including humans, form social ties and their decisions are influenced by who they are connected with and the attributes of their peers. It is not only the individual entities that are important, but the networks surrounding them hold equal significance. As we collect ever larger amounts of multi-faceted data, the opportunity and necessity to understand the world in a way that respects such interconnection becomes increasingly vital for expanding our understanding of these various fields.

In economics and social sciences, studying individual heterogeneity is an essential topic for empirical research and policy making. The inherent differences among individuals, which encompass a wide range of factors such as personal characteristics, socioeconomic backgrounds, and life circumstances, can significantly impact how individuals respond to economic policies and play a fundamental role for policy design. The literature has a long history of recognizing and accounting for individual heterogeneity, with approaches such as the demand estimation method for differentiated products (berry1995automobile), the individual fixed-effects panel-data method (wooldridge2010econometric), and the potential outcomes framework for heterogeneous treatment effects (rubin1974estimating), among others. A more recent stream of literature has started utilizing machine learning techniques to incorporate individual heterogeneity. Examples include belloni2014inference using the Lasso, bonhomme2015grouped using k-means, and wager2018estimation using random forests. In this trend, the study most closely related to our work is farrell2021deep who have demonstrated a promising potential of using artificial neural networks to capture individual heterogeneity.

The Multi-layer Perceptron (MLP) architecture, analyzed in farrell2021deep, is a potent tool to flexibly represent the dependence of individual heterogeneity on a large number of observed characteristics. However, in various empirical settings, we often encounter rich network data in conjunction with these observed characteristics. Adapting the canonical neural network model to effectively capture the variations embedded in network information across observations, which we refer to as network heterogeneity, could be a compelling topic to explore in economic research.

With this goal in mind, it is inspiring to reflect on the success of Convolutional Neural Networks (CNNs) in the deep learning literature (lecun1998gradient) and connect that to our goal of incorporating network heterogeneity. CNNs have achieved significant success in computer vision, with notable applications such as image classification, object detection, and semantic segmentation. For example, a remarkable application of CNNs is the renowned artificial intelligence program, AlphaGo, where board positions in the game of Go are represented as visual patterns and CNNs effectively detect and classify these patterns. A prominent feature in the domain of computer vision is the importance of small patterns within an image, which often matter more than the whole image in many applications. And these small patterns often appear repeatedly within and across images. The convolutional layer in CNNs consists of filters that parsimoniously extracts only the local information around a focal pixel. And the parameter sharing among these filters, meaning the same parameters are used across filters, allows for the detection of repeated patterns, contributing to the success of the CNN architecture.

Analogously, in the context of network heterogeneity, it is intuitive that peers closely connected via networks would provide more useful information than distant nodes. Also, the same aggregation function could prove beneficial for extracting local information across different nodes. Hence, it might be worthwhile to borrow insights from CNNs when studying network heterogeneity. However, a direct application is not viable, as traditional CNNs are designed to work on data residing in regular, grid-like structures, rather than on graphs.

Graph Neural Networks (GNNs) offer a general modeling framework for building neural networks on graphs, a topic thoroughly overviewed by hamilton2020graph. Intuitively, the network heterogeneity of a focal node depends on the conditions of their connected peers, which in turn depend on the conditions of their peers, and so on. GNNs provide a flexible framework for incorporating peer information, taking into account network topology. The canonical MLP is less suited to this task, as each focal node has a different number of peers with varying distances, while the MLP can only incorporate a fixed number of inputs. GNNs already have many impressive real-world applications. For instance, stokes2020deep modeled chemical molecule structures as graphs, based on which they predict their pathogen inhibitory properties using GNNs. Their study successfully identified a new molecule named halicin against Acinetobacter baumanni, which is one of the highest priority pathogens the World Health Organization urgently seeks new antibiotics for. Despite the impactful applications of GNNs in many fields, to the best of our knowledge, no study has yet employed GNNs to quantify network heterogeneity for causal inference and policy recommendations in economics. Our study aims to contribute to this area.

The primary theoretical contribution of our study lies in providing the convergence rate of the GNN estimator. The novelty of our theoretical development resides in bounding the complexity measure and approximation error, both of which are specific to GNNs, and in adapting the localization analysis (bartlett2005local) to incorporate dependent data using a dependency graph. These elements are crucial in establishing the convergence rate of the GNN estimator.

We then demonstrate the application of the GNN estimator in semiparametric causal inference. The primitive inputs for conducting robust inference of causal effects are the conditional expectations of potential outcomes and the propensity score. The GNN estimator, by effectively integrating local neighborhood information, offers flexible estimations for both the conditional expectations of potential outcomes and the propensity score, thereby facilitating causal inference. The convergence rate result of the GNN estimator elucidates its applicability in semiparametric causal inference. In broader terms, our theoretical analysis sheds light on a general procedure that future researchers could adopt to study other artificial neural network architectures tailored to their empirical needs, and subsequently conduct causal inference.

In our empirical application, we utilize the data from BanerjeeChandrasekharDufloJackson2013 to analyze individual decisions on microfinance participation given their neighborhood surroundings. Our treatment effects estimates, based on first-stage GNN estimates, indicate that well-connected households, as measured by network centrality measures, are less inclined to borrow via microfinance. Given the empirical evidence that participating households have a much higher chance to disseminate information to their neighbors than non-participating households, our findings highlight a tradeoff for information diffusion: while better-connected households have more avenues to spread information, their lower participation rates curb this potential. Using the participation probabilities predicted by the GNN estimates, we find considerable potential for improvement in microfinance information targeting by achieving a balance between the participation rate and network centrality of targeted households. This refined approach can be leveraged to facilitate information diffusion through social networks.

We organize the paper as follows. In Section (ref), we introduce the modeling setup and GNN estimator. We then present the rate of convergence for the GNN estimator in Section (ref), followed by semiparametric causal inference for treatment effects estimators in Section (ref). Monte Carlo simulations and empirical applications are provided in Sections (ref) and (ref), respectively. Finally, we conclude in Section (ref). All the theoretical proofs are contained in the appendix.

We use the following notations throughout the paper. A bold capital letter (e.g., $\mathbf{A}$) represents a matrix, a bold lowercase letter (e.g., $\mathbf{a}$) signifies a vector or vector-valued function, and an unbold letter (e.g., $A$ or $a$) denotes a scalar, unless stated otherwise. For a real number $\alpha,$ $\text{sgn}\left(\alpha\right)$ takes the value $1$ if $\alpha\geq0$ and $0$ otherwise. $\log$ refers to the natural logarithm. For vectors $\mathbf{a}$ and $\boldsymbol{b}$, $\mathbf{a}\cdot\boldsymbol{b}$ symbolizes the inner product. The letter $C$ is designated to denote a fixed finite positive constant, which does not depend on the sample size $n$. For concise presentation, the specific value of $C$ may vary across contexts, even from line to line throughout the proof. For any non-random non-negative sequences $a_{n}$ and $b_{n}$, $a_{n}=O\left(b_{n}\right)$ or $a_{n}\lesssim b_{n}$ means there exists a fixed finite constant $C$ such that $a_{n}\leq Cb_{n}$ for all $n$, and $a_{n}\asymp b_{n}$ means $a_{n}=O\left(b_{n}\right)$ and $b_{n}=O\left(a_{n}\right)$. $\mathbb{Z}_{+}$ denotes the set of positive integers. For any two real numbers $a$ and $b$, $a\vee b=\max\left\{ a,b\right\} $. For any positive integer $n$, $\left[n\right]=\left\{ 1,...,n\right\} $. And for any function $f$, $\left\Vert f\right\Vert _{\infty,\mathcal{A}}=\max_{x\in\mathcal{A}}\left|f\left(x\right)\right|$.

Graph neural networks

In this section, we first introduce the definition of network heterogeneity, and then present the GNN estimator for estimating the network heterogeneity.

Network heterogeneity

For each node $i\in\left[n\right]$, let $y_{i}\in\mathbb{R}$ denote an outcome variable and $\boldsymbol{x}_{i}\in\mathbb{R}^{d}$ a finite-dimensional vector of covariates. The adjacency matrix, denoted as $\boldsymbol{D}=\left(d_{ij}\right)$, is an $n\times n$ matrix indicating the connections between the nodes. In this paper, we focus on binary adjacency matrices, which are applicable to both undirected and directed graphs. For an undirected graph, the adjacency matrix is symmetric, where $d_{ij}=d_{ji}=1$ if there is an edge between nodes $i$ and $j$, and $d_{ij}=d_{ji}=0$ otherwise. In contrast, for a directed graph, the adjacency matrix $\boldsymbol{D}$ is normally asymmetric, where $d_{ij}=1$ if there is a directed edge from node $i$ to node $j$, and $d_{ij}=0$ if there is no such edge. We can also extend the framework to non-binary adjacency matrices to accommodate more complex information such as edge weights and labels, albeit at the expense of additional notation. Using the binary adjacency matrix, denote the set of $i$'s adjacent neighbors as $\mathcal{N}\left(i\right)=\left\{ j\in\left[n\right]:d_{ij}=1\right\} $. Also, let $\left|\mathcal{N}\left(i\right)\right|$ be the number of elements in set $\mathcal{N}\left(i\right)$, and $\left|\mathcal{N}\left(i\right)\right|=0$ if $\mathcal{N}\left(i\right)$ is empty. The researcher observes $\left\{ y_{i}\right\} _{i\in\left[n\right]}$, $\left\{ \boldsymbol{x}_{i}\right\} _{i\in\left[n\right]}$, and $\boldsymbol{D}$ for a realized sample of size $n$.

In the following, we first introduce the intuition of one-layer GNNs, and then progress to multi-layer GNNs. The type of network heterogeneity that one-layer GNNs aim to encapsulate can be expressed as

{ \[ z_{*i}\left(f_{*}\right)=f_{*}\left(\boldsymbol{x}_{i},\frac{1}{\left|\mathcal{N}\left(i\right)\right|}\sum_{j\in\mathcal{N}\left(i\right)}\boldsymbol{x}_{j}\right), \] }wherein the network heterogeneity for node $i$ depends on its own attributes, $\boldsymbol{x}_{i}$, and the average of its neighbors' attributes. The function $f_{*}$ represents an unrestricted mapping. Hence, the popular linear-in-means model, in which the outcome variable is a linear function of $\boldsymbol{x}_{i}$ and the average covariates of neighbors, can be viewed as a specific example of this representation. While our primary focus in this study is on the setup using neighbors' average, alternative aggregation methods to combine neighbors' information, such as taking the maximum, minimum, sum, or weighted sum given the availability of weights, can be easily adopted to replace the function of averaging. In this paper, we assume one of these simple aggregation functions is valid, which is a plausible assumption in many empirical contexts, and the exploration of more complex aggregation methods is deferred to future research.

While the one-layer model is straightforward, a significant limitation lies in the possible inadequacy of using the adjacent neighbors' attributes alone to summarize a node's distinctiveness. Therefore, it could be more fitting to consider a two-layer extension in which node $i$'s network heterogeneity relies on some latent embeddings of itself and its neighbors:

{ \[ z_{*i}\left(f_{*}\right)=f_{*}^{(2)}\left(\boldsymbol{h}_{*i},\frac{1}{\left|\mathcal{N}\left(i\right)\right|}\sum_{j\in\mathcal{N}\left(i\right)}\boldsymbol{h}_{*j}\right), \] }where $\boldsymbol{h}_{*i}\in\mathbb{R}^{d_{h*}}$ represent some unobserved measures of node $i$'s uniqueness with an unknown dimension $d_{h*}\in\mathbb{Z}_{+}$. This concept is depicted in Figure (ref). The hidden embeddings of node $i$ consequently depend on the observed covariates of both itself and its neighbors:

{ \[ \boldsymbol{h}_{*i}=\boldsymbol{f}_{*}^{(1)}\left(\boldsymbol{x}_{i},\frac{1}{\left|\mathcal{N}\left(i\right)\right|}\sum_{j\in\mathcal{N}\left(i\right)}\boldsymbol{x}_{j}\right). \] }

In this two-layer model, the network heterogeneity of node $i$, denoted as $z_{*i}\left(f_{*}\right)$, depends on its own attributes, its neighbors' attributes, and the attributes of its neighbors' neighbors, utilizing the local network structure. This is a natural extension of the one-layer model setup. In many empirical contexts, an individual's decision (for example, taking a loan) depends not only on the economic status of her immediate connections, as their willingness to lend may also depend on who they are connected with and could borrow from. This concept can intuitively be extended to setups with multiple layers. For instance, in a three-layer model, the third-degree connections (i.e., neighbors of neighbors' neighbors) can also be relevant in an individual's decision-making process.

To elucidate the idea more formally, we define an $L$-layer GNNs model iteratively as $\forall l\in\left[L\right]$,{ \[ \boldsymbol{h}_{*i}^{\left(l\right)}\left(\boldsymbol{f}\right)=\boldsymbol{f}^{\left(l\right)}\left(\boldsymbol{h}_{*i}^{\left(l-1\right)}\left(\boldsymbol{f}\right),\overline{\boldsymbol{h}}_{*\mathcal{N}\left(i\right)}^{\left(l-1\right)}\left(\boldsymbol{f}\right)\right)\in\mathbb{R}^{d_{h*}^{\left(l\right)}}, \] }where $\boldsymbol{h}_{*i}^{\left(0\right)}\left(\boldsymbol{f}\right)=\boldsymbol{x}_{i}$, $\overline{\boldsymbol{h}}_{*\mathcal{N}\left(i\right)}^{\left(l-1\right)}\left(\boldsymbol{f}\right)=\frac{1}{\left|\mathcal{N}\left(i\right)\right|}\sum_{j\in\mathcal{N}\left(i\right)}\boldsymbol{h}_{*j}^{\left(l-1\right)}\left(\boldsymbol{f}\right)$, and $\overline{\boldsymbol{h}}_{*\mathcal{N}\left(i\right)}^{\left(l-1\right)}\left(\boldsymbol{f}\right)=\boldsymbol{0}$ if $\mathcal{N}\left(i\right)$ is empty, with $\boldsymbol{f}$ being the set of functions constructing the final layer output, i.e., $\boldsymbol{f}=\left\{ \boldsymbol{f}^{\left(l\right)}\right\} {}_{l\in\left[L\right]}$ and $\boldsymbol{f}^{\left(l\right)}=\left\{ f_{j}^{\left(l\right)}\right\} {}_{j\in\left[d_{h*}^{\left(l\right)}\right]}$. The latent embedding dimension in each intermediate layer $l\in\left[L-1\right]$ is an unknown finite fixed integer $d_{h*}^{\left(l\right)}\in\mathbb{Z}_{+}$, and the final layer output is a scalar ($d_{h*}^{\left(L\right)}=1$). We denote the final layer output of the $L$-layer GNNs model as $z_{*i}\left(\boldsymbol{f}\right)\coloneqq\boldsymbol{h}_{*i}^{\left(L\right)}\left(\boldsymbol{f}\right)\in\mathbb{R}$. Hence, the formulation of $z_{*i}\left(\boldsymbol{f}\right)$ depends on the number of layers, $L$, and the dimensions of latent embeddings in each intermediate layer, $\left\{ d_{h*}^{\left(l\right)}\right\} _{l\in\left[L-1\right]}$. To simplify notation, we adopt the convention in this paper to omit the reference to $L$ and $\left\{ d_{h*}^{\left(l\right)}\right\} _{l\in\left[L-1\right]}$ in $z_{*i}\left(\boldsymbol{f}\right)$. Also, define $\boldsymbol{f}_{*}$ as a set of functions yielding the smallest average loss function value,

{

}Common examples of the loss function $\ell$ include the least squares loss $\ell\left(y,z\right)=\frac{1}{2}\left(y-z\right)^{2}$ for real-valued outcome variables, and the negative log-likelihood of logistic regression $\ell\left(y,z\right)=-yz+\log\left(1+\exp\left(z\right)\right)$ for binary outcome variables.

We call the final layer output associated with an optimal $\boldsymbol{f}_{*}$, namely $z_{*i}\left(\boldsymbol{f}_{*}\right)$, the network heterogeneity of node $i$. On average, the network heterogeneity yields the smallest loss function value for predicting the outcome variable $y_{i}$ given the local neighborhood information. Standard derivations show that $z_{*i}\left(\boldsymbol{f}_{*}\right)=\mathbb{E}\left[y_{i}\mid\boldsymbol{h}_{*i}^{\left(L-1\right)}\left(\boldsymbol{f}_{*}\right),\overline{\boldsymbol{h}}_{*\mathcal{N}\left(i\right)}^{\left(L-1\right)}\left(\boldsymbol{f}_{*}\right)\right]$ for the least squares loss, and $\frac{\exp\left(z_{*i}\left(\boldsymbol{f}_{*}\right)\right)}{1+\exp\left(z_{*i}\left(\boldsymbol{f}_{*}\right)\right)}=\mathbb{E}\left[y_{i}\mid\boldsymbol{h}_{*i}^{\left(L-1\right)}\left(\boldsymbol{f}_{*}\right),\overline{\boldsymbol{h}}_{*\mathcal{N}\left(i\right)}^{\left(L-1\right)}\left(\boldsymbol{f}_{*}\right)\right]$ for the logistic loss. In other words, the network heterogeneity determines the conditional expectation of the outcome variable given the local neighborhood information as summarized in the penultimate layer latent embeddings. Without further restrictions, $z_{*i}\left(\boldsymbol{f}_{*}\right)$ may not be unique, while our main theoretical result in Theorem (ref) is valid for any $z_{*i}\left(\boldsymbol{f}_{*}\right)$ that corresponds to such an optimal forecasting rule.

The concept of network heterogeneity, represented as $z_{*i}\left(\boldsymbol{f}_{*}\right)$, holds potential for a wide range of applications in empirical research. This paper focuses on a significant application area: the semiparametric causal inference of treatment effects. Central to identifying and robustly inferring causal effects are the conditional expectations of potential outcomes and the propensity score, both could depend on local network surroundings. As detailed in Section (ref), under appropriate assumptions, the conditional expectation of the potential outcome $y_{i}\left(t\right)$ and the propensity score are known functions of respective network heterogeneity variables, $z_{*i}\left(\boldsymbol{f}_{*}^{t}\right)$ and $z_{*i}\left(\boldsymbol{f}_{*}^{p}\right)$. These network heterogeneities enable flexible characterization of dependencies on local network surroundings. Utilizing GNN estimators, we effectively estimate these network heterogeneities, leading to accurate estimates for both the conditional expectations of potential outcomes and the propensity score. These estimates are crucial for conducting robust inference on various causal effects. In Section (ref), we focus on a specific causal effect parameter: the average effect of a counterfactual policy $\pi\left(s\right)=\mathbb{E}\left[s\left(\boldsymbol{\xi}_{i}\right)y_{i}\left(1\right)+\left(1-s\left(\boldsymbol{\xi}_{i}\right)\right)y_{i}\left(0\right)\right]$, which plays an integral role in our empirical application.

Beyond the estimation of treatment effects, another noteworthy application area of the concept of network heterogeneity, $z_{*i}\left(\boldsymbol{f}_{*}\right)$, lies in improving structural estimations. Building upon and extending the research by farrell2021individual, the observation-specific coefficients in structural models can be formulated as functions of the observation-specific local network data via the concept of network heterogeneity. This flexibility allows, for instance, the individual-specific consumer demand elasticities to vary over consumers' social network scenarios in demand estimations. As a result, the network heterogeneity enriched structural models can be powerful to study key economic parameters, such as elasticity and surplus, and to answer policy questions, such as optimal pricing and targeting. The empirical exploration of these diverse applications presents an opportunity for future research.

We now introduce some additional notations to describe network topology. First, it would be convenient to transform the set $\mathcal{N}\left(i\right)$ into a tuple (i.e., an ordered list which may contain multiple occurrences of the same element). Without loss of generality, we could fix the labels of nodes in a sample using natural numbers. Then, let $\mathcal{D}\left(i\right)$ denote the tuple of nodes in the set $\mathcal{N}\left(i\right)$ arranged according to the order of natural numbers. For instance, set $\mathcal{D}\left(i\right)=\left(1,2,3\right)$ if $\mathcal{N}\left(i\right)=\left\{ 1,2,3\right\} $ and $\mathcal{D}\left(i\right)=\emptyset$ if $\mathcal{N}\left(i\right)=\emptyset$. Also, denote the concatenation of two tuples as $\left(a_{1},...,a_{m}\right)\oplus\left(b_{1},...,b_{n}\right)=\left(a_{1},...,a_{m},b_{1},...,b_{n}\right)$.\footnote{If there are empty tuples, set $\left(a_{1},...,a_{m}\right)\oplus\emptyset=\left(a_{1},...,a_{m}\right)$, $\emptyset\oplus\left(b_{1},...,b_{n}\right)=\left(b_{1},...,b_{n}\right)$, $\emptyset\oplus\emptyset=\emptyset$, and $\oplus_{j\in\emptyset}\mathcal{D}\left(j\right)=\emptyset$.} Then, let $\mathcal{D}_{l}\left(i\right)$ denote the tuple of nodes whose distance from node $i$ is $l$ for $l\in\left\{ 0,1,...,L\right\} $. In particular, define $\mathcal{D}_{l}\left(i\right)$ recursively such that $\mathcal{D}_{0}\left(i\right)=\left(i\right)$ and $\mathcal{D}_{l}\left(i\right)=\oplus_{j\in\mathcal{\mathcal{D}}_{l-1}\left(i\right)}\mathcal{D}\left(j\right)$. Using the example in Figure (ref) as an illustration, we have that $\mathcal{D}_{0}\left(i\right)=\left(i\right)$, which includes node $i$ itself; $\mathcal{D}_{1}\left(i\right)=\mathcal{D}\left(i\right)=\left(1,2,3\right)$, containing the immediate neighbors of node $i$; and $\mathcal{D}_{2}\left(i\right)=\oplus_{j\in\left(1,2,3\right)}\mathcal{D}\left(j\right)=\left(4,5,6,7,8,9\right)$, including the neighbors of node $i$'s neighbors.

Moreover, let $\boldsymbol{\xi}_{i,l}$ denote the local network information for node $i$ up to distance $l$ for $l\in\left\{ 0,1,...,L\right\} $. Specifically, define $\boldsymbol{\xi}_{i,l}$ recursively such that $\boldsymbol{\xi}_{i,0}=\left(\boldsymbol{x}_{i}\right)$ and $\boldsymbol{\xi}_{i,l}=\left(\boldsymbol{\xi}_{i,l-1},\left(\left(\boldsymbol{x}_{k}\right)_{k\in\mathcal{D}\left(j\right)}\right)_{j\in\mathcal{D}_{l-1}\left(i\right)}\right)$.\footnote{When encountering the empty tuple $\emptyset$, set $\left(\boldsymbol{x}_{k}\right)_{k\in\emptyset}=\emptyset$ and $\left(\left(\boldsymbol{x}_{k}\right)_{k\in\mathcal{D}\left(j\right)}\right)_{j\in\emptyset}=\left(\emptyset\right)$.} In the example in Figure (ref), it is easy to find that $\boldsymbol{\xi}_{i,0}=\left(\boldsymbol{x}_{i}\right)$, $\boldsymbol{\xi}_{i,1}=\left(\boldsymbol{\xi}_{i,0},\left(\left(\boldsymbol{x}_{1},\boldsymbol{x}_{2},\boldsymbol{x}_{3}\right)\right)\right),$ and $\boldsymbol{\xi}_{i,2}=\left(\boldsymbol{\xi}_{i,1},\left(\left(\boldsymbol{x}_{4},\boldsymbol{x}_{5}\right),\left(\boldsymbol{x}_{6},\boldsymbol{x}_{7}\right),\left(\boldsymbol{x}_{8},\boldsymbol{x}_{9}\right)\right)\right).$ This definition of $\boldsymbol{\xi}_{i,l}$ can accommodate scenarios where nodes share common friends and networks have cycles (see Figure (ref) for illustrations). For simplicity in notation, we denote the $L$-hop network information from the perspective of node $i$ as $\boldsymbol{\xi}_{i}\coloneqq\boldsymbol{\xi}_{i,L}$ (omitting reference to $L$). With this definition of $\boldsymbol{\xi}_{i}$, it is evident that $\boldsymbol{\xi}_{i}$ encapsulates all the needed input data to formulate the network heterogeneity variable $z_{*i}\left(\boldsymbol{f}_{*}\right)$ for node $i$.

GNN estimator

We now introduce the GNN estimator, which is utilized to estimate the network heterogeneity variable $z_{*i}\left(\boldsymbol{f}_{*}\right)$. In particular, for each layer $l\in\left[L\right]$, the latent feature embeddings are updated recursively using shallow neural networks as { \[ \boldsymbol{h}_{i}^{\left(l\right)}=\mathbf{\bm{\sigma}}\left(\boldsymbol{A}^{\left(l\right)}\boldsymbol{h}_{i}^{\left(l-1\right)}+\boldsymbol{A}_{\mathcal{N}}^{\left(l\right)}\overline{\boldsymbol{h}}_{\mathcal{N}\left(i\right)}^{\left(l-1\right)}+\boldsymbol{b}^{\left(l\right)}\right)\in\mathbb{R}^{d_{h}^{\left(l\right)}}, \] }where $\bm{\sigma}$ is an activation function applied element-wise, and $\left(\boldsymbol{A}^{\left(l\right)},\boldsymbol{A}_{\mathcal{N}}^{\left(l\right)},\boldsymbol{b}^{\left(l\right)}\right)$ are parameters to be estimated. The initial feature embeddings are defined by $\boldsymbol{h}_{i}^{(0)}=\boldsymbol{x}_{i}$, and the mean feature embeddings of neighboring nodes are $\overline{\boldsymbol{h}}_{\mathcal{N}\left(i\right)}^{\left(l-1\right)}=\frac{1}{\left|\mathcal{N}\left(i\right)\right|}\sum_{j\in\mathcal{N}\left(i\right)}\boldsymbol{h}_{j}^{\left(l-1\right)}$ (if $\mathcal{N}\left(i\right)$ is empty, set $\overline{\boldsymbol{h}}_{\mathcal{N}\left(i\right)}^{\left(l-1\right)}=\boldsymbol{0}$). Then, the feature of node $i$ is obtained through the transformation

{

}where $\left(\boldsymbol{a},b\right)$ are additional parameters, and $\boldsymbol{\theta}$ is the vectorization of all parameters appeared in the construction of $z_{i}\left(\boldsymbol{\theta}\right)$. Note that, given $\boldsymbol{\theta}$, $z_{i}\left(\boldsymbol{\theta}\right)$ can be fully constructed using the input data $\boldsymbol{\xi}_{i}$.

Denote the parameter space of $\boldsymbol{\theta}$ by $\Theta_{d_{h}}$, which is indexed by the dimensions of the embeddings for each layer, i.e., $d_{h}=\left\{ d_{h}^{\left(l\right)}\right\} _{l\in\left[L\right]}$. In particular,{ \[ \Theta_{d_{h}}=\Theta_{\left\{ d_{h}^{\left(1\right)},...,d_{h}^{\left(L\right)}\right\} }=\left\{ \boldsymbol{a}\in\mathbb{R}^{d_{h}^{\left(L\right)}},b\in\mathbb{R},\left\{ \boldsymbol{A}^{\left(l\right)},\boldsymbol{A}_{\mathcal{N}}^{\left(l\right)}\in\mathbb{R}^{d_{h}^{\left(l\right)}\times d_{h}^{\left(l-1\right)}},\boldsymbol{b}^{\left(l\right)}\in\mathbb{R}^{d_{h}^{\left(l\right)}}\right\} _{l\in\left[L\right]}\right\} , \] }with $d_{h}^{\left(0\right)}=d$ being the number of covariates.

Moreover, denote $\left\Vert z_{i}\left(\boldsymbol{\theta}\right)\right\Vert _{\infty}=\underset{\boldsymbol{\xi}_{i}}{\text{sup}}\left|z_{i}\left(\boldsymbol{\theta}\right)\right|$, which, given $\boldsymbol{\theta}$, is the largest value of $\left|z_{i}\left(\boldsymbol{\theta}\right)\right|$ over all possible $L$-hop local networks $\boldsymbol{\xi}_{i}$. The value of $\left\Vert z_{i}\left(\boldsymbol{\theta}\right)\right\Vert _{\infty}$ depends on $\boldsymbol{\theta}$ and is common across all nodes. We further refine the set $\Theta_{d_{h}}$ so that $\left\Vert z_{i}\left(\boldsymbol{\theta}\right)\right\Vert _{\infty}$ is bounded from above, which results in a new set defined as{ \[ \Theta_{d_{h},\bar{z}}=\left\{ \boldsymbol{\theta}\in\Theta_{d_{h}}:\left\Vert z_{i}\left(\boldsymbol{\theta}\right)\right\Vert _{\infty}\leq\bar{z}\right\} , \] }where $\bar{z}$ is a fixed constant used throughout the analysis.

To estimate the parameters $\boldsymbol{\theta}$, we consider the optimization problem{

}with $\ell$ being the same loss function used in the formulation of $z_{*i}\left(\boldsymbol{f}_{*}\right)$ in ((ref)). We refer to $z_{i}\left(\hat{\boldsymbol{\theta}}\right)$ as the GNN estimator of the network heterogeneity $z_{*i}\left(\boldsymbol{f}_{*}\right)$. The minimizer $\hat{\boldsymbol{\theta}}$ and hence $z_{i}\left(\hat{\boldsymbol{\theta}}\right)$ may not be unique, while our main theoretical result in Theorem (ref) is valid for any minimizer $\hat{\boldsymbol{\theta}}$ that solves the problem in ((ref)). Without loss of generality, we let $\bar{z}>1.1M$ with $M$ being the upper bound of $\left|z_{*i}\left(\boldsymbol{f}_{*}\right)\right|$ as imposed in Assumption I 1 below (e.g., we could set $\bar{z}=2M$ as in farrell2021deep). In the following, we present the theoretical properties of the GNN estimator $z_{i}\left(\hat{\boldsymbol{\theta}}\right)$ to estimate the target object $z_{*i}\left(\boldsymbol{f}_{*}\right)$.

Theoretical property

In this section, we present the convergence rate of the GNN estimator. We first define the Sobolev space { \[ \mathcal{W}_{\eta}^{\beta,\infty}\left(\Omega\right)=\left\{ f:\sum_{0\leq|\boldsymbol{k}|\leq\beta}\underset{

}{\ensuremath{ess\ sup}}\left|\frac{\partial^{|\mathbf{k}|}f}{\partial x_{1}^{k_{1}}...\partial x_{p}^{k_{p}}}\right|\leq\eta\right\} , \] }where $\Omega\subseteq\mathbb{R}^{p}$, $\boldsymbol{k}\coloneqq\left(k_{1},...,k_{p}\right)$ with $k_{1},...,k_{p}$ being nonnegative integers, $\left|\boldsymbol{k}\right|=k_{1}+...+k_{p}$, and $\boldsymbol{x}=\left(x_{1},...,x_{p}\right)$ being an argument of $f$. Also, we define{ \[ \boldsymbol{\theta}_{*}\in\underset{

}{\ensuremath{arg\min}}\mathbb{E}\left[\ell\left(y_{i},z_{i}\left(\boldsymbol{\theta}\right)\right)\right]. \] }

We now introduce the assumptions to demonstrate the convergence rate.

Assumption I

Assumption I 1 postulates that the functions composing network heterogeneity are smooth with bounded input and output spaces. Given the bounded nature of these input and output spaces, it is innocuous to normalize the input and output spaces of each internal layer to be unit cubes. Our theory can be adapted to allow the dimensions of hidden embeddings, $d_{h*}^{\left(1\right)},...,d_{h*}^{\left(L-1\right)}$, to grow slowly with the sample size, but we omit that complexity for concise presentation.

Assumption I 2 accommodates common activation functions such as the sigmoid, $\sigma\left(x\right)=\frac{1}{1+\exp\left(-x\right)}$, which takes four operations, and the tanh, $\sigma\left(x\right)=\frac{\exp\left(x\right)-\exp\left(-x\right)}{\exp\left(x\right)+\exp\left(-x\right)}$, which takes six operations. As this paper focuses on the most prevalently used GNN architectures, where each layer comprises a shallow neural network, our setup does not account for activation functions that are particularly effective for deep neural networks, such as the ReLU, which is non-differentiable at zero.

Assumption I 3 posits that the curvature of the loss function is bounded from below and above at $z_{*i}\left(\boldsymbol{f}_{*}\right)$, and bounded from below at $z_{i}\left(\boldsymbol{\theta}_{*}\right)$. Furthermore, the loss function is Lipschitz continuous with respect to the second argument. Regarding ((ref)) and ((ref)), suppose $\mathbb{E}\left[y_{i}\mid\boldsymbol{h}_{*i}^{\left(L-1\right)}\left(\boldsymbol{f}_{*}\right),\overline{\boldsymbol{h}}_{*\mathcal{N}\left(i\right)}^{\left(L-1\right)}\left(\boldsymbol{f}_{*}\right)\right]=\mathbb{E}\left[y_{i}\mid\boldsymbol{\xi}_{i}\right]$, which indicate that the penultimate layer latent embeddings sufficiently summarize the local neighborhood information needed for the conditional mean of the outcome variable. Then, farrell2021deep shows that ((ref)) and ((ref)) hold true for many commonly used loss functions, including the least squares loss for bounded real-valued outcome variables, and logistic loss for binary outcome variables. Condition ((ref)) is introduced in Section 5.2 of bartlett2005local with valid examples. If we could show that $\hat{\boldsymbol{\theta}}$ is a consistent estimator of $\boldsymbol{\theta}_{*}$, we can relax this condition by assuming it holds only for $\boldsymbol{\theta}\in\Theta_{d_{h},\bar{z}}$ in a small neighborhood around $\boldsymbol{\theta}_{*}$, and this weaker condition can be easily verified for any loss function under mild regularity conditions. Meanwhile, the consistency of the M-estimator $\hat{\boldsymbol{\theta}}$ can be established using Theorem 5.7 in \Citet{Vaart1998}, where the uniform convergence of the criterion function can be demonstrated. Yet, to simply the proof, we employ the stronger condition in this paper as in bartlett2005local.

Under Assumption I 4, $\left(y_{i},\boldsymbol{\xi}_{i}\right)$ is identically distributed over $i\in\left[n\right]$, and each node in the sample of size $n$ has at most $c_{n}$ neighbors. It would be more precise to use the triangular array notation such that for each $n$, $\left\{ \left(y_{n,i},\boldsymbol{\xi}_{n,i}\right)\right\} _{i\in\left[n\right]}$ denotes a sequence of identically distributed set of random variables each distributed as $\left(y_{n,i},\boldsymbol{\xi}_{n,i}\right)\sim P_{n}$. To streamline our notation, we suspend the use of subscript $n$ in $\left(y_{n,i},\boldsymbol{\xi}_{n,i}\right)$ and in any object depending on them ($c_{n}$ is an exception), although our analysis is valid for the setting where the distribution of $\left(y_{n,i},\boldsymbol{\xi}_{n,i}\right)$ varies with $n$.

Assumption I 4 does not require the set of random variables $\left(y_{i},\boldsymbol{\xi}_{i}\right)$ to be independent across $i$. In our context, accounting for dependence is important. Examples of the sources of dependence include (1) node attributes $\boldsymbol{x}_{i}$ could be dependent across $i$, (2) the $L$-hop local neighborhood $\boldsymbol{\xi}_{i}$ may overlap across $i$, and (3) the edge formation process itself might induce dependence. To account for dependence, we utilize the concepts of dependency graph and smallest proper cover, as introduced below.

Given the dependency graph of $\left\{ \left(y_{i},\boldsymbol{\xi}_{i}\right)\right\} _{i\in\left[n\right]}$, we could construct the smallest proper cover based on the dependency graph denoted by $\left\{ \mathcal{C}_{j}\right\} _{j\in\left[J\right]}$.\footnote{If the smallest proper cover is not unique, selecting any one of them as $\left\{ \mathcal{C}_{j}\right\} _{j\in\left[J\right]}$ would suffice.} By construction, $\left(y_{i},\boldsymbol{\xi}_{i}\right)$'s are independent within each cover $\mathcal{C}_{j}$ but can be arbitrarily correlated across covers. Both the number of covers $J$ and the set of covers $\left\{ \mathcal{C}_{j}\right\} _{j\in\left[J\right]}$ depend on the sample size $n$, although we leave out the subscript $n$ to simplify notations. In our analysis, we treat $\left\{ \mathcal{C}_{j}\right\} _{j\in\left[J\right]}$ as fixed for each sample size $n$, otherwise the analysis becomes overly complicated.

The construction of the smallest proper cover based on dependency graph can be better understood with the help of graph coloring, as discussed in ralaivola2015entropy. The idea of graph coloring is to assign different colors to nodes in a graph so that no two adjacent nodes share the same color. Using the dependency graph for graph coloring, nodes assigned to the same color are necessarily independent, thereby enabling the construction of a proper cover in which each cover $\mathcal{C}_{j}$ includes only the nodes of the same color $j$. Note that this graph coloring is based on the dependency graph, rather than the observed network. Subsequently, a smallest proper cover can be obtained through graph coloring using the minimal number of colors needed.\footnote{While the concept of graph coloring helps to understand the theoretical results, in practice there is no necessity for empirical researchers to construct the smallest proper cover for inference, as demonstrated in Section (ref) for the semiparametric inference.}

Remark {[}Motivating example{]}: To illustrate, consider the setting of our empirical application, where the dataset includes many villages. In this setting, we suppose that households within the same village can form links and have correlated attributes and outcome variables. However, households across different villages are unable to form links and their attributes and outcome variables are considered independent. This leads to a specific pattern in our dependency graph for the set $\left\{ \left(y_{i},\boldsymbol{\xi}_{i}\right)\right\} _{i\in\left[n\right]}$: households within each village are interconnected, forming a complete graph, while there are no connections between households from different villages.

Given that households in the same village are completely connected in the dependency graph, the graph coloring principle dictates that each household in the same village must be assigned a unique color. Hence, the number of covers of a smallest proper cover ($J$), which equates to the smallest number of colors needed, corresponds to the size of the largest village. Moreover, the largest cover size $\left(\max_{j\in\left[J\right]}\left|\mathcal{C}_{j}\right|\right)$ cannot exceed the number of villages. With appropriate coloring, we could set $\left|\mathcal{C}_{j}\right|\asymp\frac{n}{J}$ for every $j\in\left[J\right]$, ensuring that each cover is similar in size.

As discussed previously, our study treats the smallest proper cover $\left\{ \mathcal{C}_{j}\right\} _{j\in\left[J\right]}$ as fixed for each sample size $n$. Additionally, we assume that $\left(y_{i},\boldsymbol{\xi}_{i}\right)$ is identically distributed over $i\in\left[n\right]$, and the maximum degree in $\boldsymbol{D}$ is bounded by $c_{n}$. Now we delve into additional details to justify this modeling choice.

First, suppose that nature assigns each household to a village. Both the number of villages and the specific assignment of households to these villages are considered fixed for each sample size, conditional on which we form our statistical analysis. Then, given the village assignment, households randomly draw their vector of observed covariates $\boldsymbol{x}_{i}$ and unobserved attributes $\boldsymbol{\epsilon}_{i}$ such that (1) $\left(\boldsymbol{x}_{i},\boldsymbol{\epsilon}_{i}\right)$'s are identically distributed across $i\in\left[n\right]$, (2) $\left(\boldsymbol{x}_{i},\boldsymbol{\epsilon}_{i}\right)$'s are independent for households in different villages, and (3) $\left(\boldsymbol{x}_{i},\boldsymbol{\epsilon}_{i}\right)$'s are exchangeable over $i$ within the same village. Next, households in the same village form network connections. This process could follow various commonly used network formation models. Examples include the Erd\H{o}s-R\'enyi model, in which each pair of nodes is connected with a constant probability, the Barab\'asi-Albert model, which introduces a preferential attachment process, and a generalized Erd\H{o}s-R\'enyi model, where the probability of edge formation between two nodes $i$ and $j$ depends on their respective attributes $\left(\boldsymbol{x}_{i},\boldsymbol{\epsilon}_{i}\right)$ and $\left(\boldsymbol{x}_{j},\boldsymbol{\epsilon}_{j}\right)$, and so on. An additional modification to the network formation model can be applied to fulfill the requirement that each household has at most $c_{n}$ connections, as long as this adjustment is identical across $i\in\left[n\right]$. Finally, the outcome variable $y_{i}$ is formed based on $\boldsymbol{\xi}_{i}$, and possibly on the unobserved attributes $\boldsymbol{\epsilon}_{i}$.

Following this procedure, if we are under the scenario where the sizes of villages are either identical or approximately similar (to the extent that the discrepancy in village size is negligible when $n$ is large), the assumption that $\left(y_{i},\boldsymbol{\xi}_{i}\right)$ is identically distributed over $i\in\left[n\right]$ holds true, or at least serves as a reasonable approximation. We defer the relaxation of the identical-distribution assumption to future research. $\square$

Rate of convergence

Utilizing the concept of the smallest proper cover $\left\{ \mathcal{C}_{j}\right\} _{j\in\left[J\right]}$, we now present the population and empirical $L_{2}$ convergence rates of the GNN estimator in Theorem (ref). The proof of this theorem can be found in Appendix (ref).

Theorem (ref) implies that the GNN estimator is consistent under relatively mild conditions that $J$ and $c_{n}$ do not grow too fast with $n$. In particular, if the growth rates of $J$ and $c_{n}$ satisfy $J\log J=o\left(n\right)$ and $c_{n}=o\left(\left(\frac{n/J}{\log\left(n/J\right)}\right)^{1/s}\right)$, then by Jensen's inequality the GNN estimator is consistent in terms of the population and empirical $L_{2}$ norms.

Theorem (ref) also sheds light on the usage of the GNN estimator for semiparametric inference. As will be shown in Section (ref), the GNN estimator can be used to conduct inference on treatment effects if its empirical $L_{2}$ convergence rate is $o_{p}(n^{-\frac{1}{4}})$. To achieve this rate, Theorem (ref) ensures that it is sufficient if $J$ and $c_{n}$ do not grow too fast with $n$ such that $J\log J=o\left(n^{1/2}\right)$ and {

}Suppose that the growth rates of $J$ and $c_{n}$ satisfy $J=O\left(n^{\alpha}\right)$ and $c_{n}=O\left(n^{\gamma}\right)$. Then Condition ((ref)) is satisfied if $\frac{kd_{h*}}{\beta}<0.999-2\alpha-2\gamma s$, where $k$ and $s$ are specified as in Theorem (ref). This condition is achievable if the growth rates of $J$ and $c_{n}$, and the dimensions of the input variables and latent embeddings (as captured by $d_{h*}$) are not too large, and the target function is sufficiently smooth.

It is worthwhile to compare the convergence rate of the GNN with that of the MLP in farrell2021deep. Given that the MLP approach in farrell2021deep is not designed for dependent data in networks, it is appropriate to set $J=O\left(1\right)$ and $c_{n}=O\left(1\right)$ for comparison. As our setting allows observations to be arbitrarily correlated across covers, setting $J=O\left(1\right)$ ensures that the effective number of independent observations in our setting is proportional to the sample size $n$, rendering the GNN and MLP setups comparable. Also, let $k=2$ to consider the GNN model with one layer, where the activation function does not involve taking exponentials, mirroring the setting in farrell2021deep, and let $d_{h*}=d$ for simplicity. Then, Condition ((ref)) is satisfied if $\frac{2}{1-\delta}d<\beta$, where $\delta$ can be any arbitrarily small constant. As a comparison, for the MLP with $2d$ input variables, which include an observation's own characteristics and the average characteristics of adjacent nodes, Theorem 3 in farrell2021deep requires that $\beta>2d$ for semiparametric inference. Hence, the requirements on the smoothness parameter $\beta$ for the GNN and MLP frameworks are comparable. Also, with more layers ($L\geq2$) and other settings the same, the smoothness condition becomes $\frac{4}{1-\delta}d<\beta$, which does not further change with $L$. This suggests that although the GNN approach could incorporate the characteristics of neighbors over multiple distances, it does not necessarily lead to a curse of dimensionality issue with larger $L$ when $L\geq2$.

In our empirical setting with many villages, let $n_{v}$ denote the number of observations in village $v$, and $\bar{v}$ the number of villages. Suppose that the village sizes are uniformly proportional to $J$ (i.e., $\min_{v\in\left[\bar{v}\right]}n_{v}\asymp J=\max_{v\in\left[\bar{v}\right]}n_{v}$), then Condition ((ref)) is equivalent to $\frac{\log\left(\bar{v}\right)}{\bar{v}}\left(c_{n}\right)^{s}=o\left(n^{-\frac{\beta+kd_{h*}}{2\beta}}\right)$. If the growth rates of $\bar{v}$ and $c_{n}$ satisfy $n^{\eta}\lesssim\bar{v}$ and $c_{n}=O\left(n^{\gamma}\right)$, then Condition ((ref)) holds when $\frac{kd_{h*}}{\beta}<2\eta-2\gamma s-1.001$, which is achievable if the number of villages grows sufficiently fast, both the growth rate of $c_{n}$ and the dimensions of the input variables and latent embeddings are not too large, and the target function is sufficiently smooth.

Semiparametric inference

In this section, we show that the theoretical guarantee of the GNN estimator in Theorem (ref) can be applied to conduct causal inference on parameters with heterogeneous treatment effects. This section closely follows farrell2021deep and hence our exposition is brief. More detailed explanations of the setup can be found in farrell2021deep and references therein. The main distinction in this section is our need to adapt the inference procedure to accommodate dependent data with networks.

Setup

Consider potential outcomes $\left(y_{i}\left(1\right),y_{i}\left(0\right)\right)$ for a binary treatment $t_{i}\in\left\{ 0,1\right\} $. Then the observable outcome can be expressed as $y_{i}=t_{i}y_{i}\left(1\right)+\left(1-t_{i}\right)y_{i}\left(0\right)$.

To illustrate the idea, we follow farrell2021deep and focus on the average outcome for a counterfactual treatment policy $\pi\left(s\right)=\mathbb{E}\left[s\left(\boldsymbol{\xi}_{i}\right)y_{i}\left(1\right)+\left(1-s\left(\boldsymbol{\xi}_{i}\right)\right)y_{i}\left(0\right)\right]$, where $s\left(\boldsymbol{\xi}_{i}\right)\in\left\{ 0,1\right\} $ is a known deterministic counterfactual policy that assigns treatment status to each observation based on $\boldsymbol{\xi}_{i}$. For convenience, we denote $s_{1}\left(\boldsymbol{\xi}_{i}\right)=s\left(\boldsymbol{\xi}_{i}\right)$ and $s_{0}\left(\boldsymbol{\xi}_{i}\right)=1-s\left(\boldsymbol{\xi}_{i}\right)$. Similar analysis directly applies to other popular parameters of interest, such as the average treatment effect (ATE) $\tau=\mathbb{E}\left[y_{i}\left(1\right)-y_{i}\left(0\right)\right]$ and the average treatment effect on the treated (ATT) $\mathbb{E}\left[y_{i}\left(1\right)-y_{i}\left(0\right)\mid t_{i}=1\right]$.

For $t\in\left\{ 0,1\right\} $, define

{

}where the loss function is the least squares loss for real-valued dependent variables and logistic loss for binary dependent variables. Then, $z_{*i}\left(\boldsymbol{f}_{*}^{t}\right)$ and $z_{*i}\left(\boldsymbol{f}_{*}^{p}\right)$ are the individual network heterogeneities of node $i$ that affect the potential outcomes and treatment assignment, respectively. The construction of the network heterogeneities $z_{*i}\left(\boldsymbol{f}_{*}^{t}\right)$ and $z_{*i}\left(\boldsymbol{f}_{*}^{p}\right)$ allows their dependence on node $i$'s local neighborhood surroundings. The optimization problems in ((ref)) are specific instances of the one in ((ref)), where we replace the generic outcome variable by the potential outcomes and treatment variable.

Then, the optimization problems in ((ref)) imply that

and

where, as usual, $\boldsymbol{h}_{*i,t}^{\left(L-1\right)}\left(\boldsymbol{f}_{*}^{t}\right)$ and $\boldsymbol{h}_{*i,p}^{\left(L-1\right)}\left(\boldsymbol{f}_{*}^{p}\right)$ are the penultimate layer latent embeddings in the construction of $z_{*i}\left(\boldsymbol{f}_{*}^{t}\right)$ and $z_{*i}\left(\boldsymbol{f}_{*}^{p}\right)$, respectively, and $\overline{\boldsymbol{h}}_{*\mathcal{N}\left(i\right),t}^{\left(L-1\right)}\left(\boldsymbol{f}_{*}^{t}\right)$ and $\overline{\boldsymbol{h}}_{*\mathcal{N}\left(i\right),p}^{\left(L-1\right)}\left(\boldsymbol{f}_{*}^{p}\right)$ are the corresponding average latent embeddings of the adjacent neighbors.

For convenience, we denote the conditional expectation of the potential outcome $y_{i}\left(t\right)$ and of the treatment status $t_{i}$ as $\mu_{t}\left(\boldsymbol{\xi}_{i}\right)$ and $p\left(\boldsymbol{\xi}_{i}\right),$ as indicated in ((ref)) and ((ref)). As such, $\mu_{t}\left(\boldsymbol{\xi}_{i}\right)$ and $p\left(\boldsymbol{\xi}_{i}\right)$ can be expressed by the network heterogeneity variables, $z_{*i}\left(\boldsymbol{f}_{*}^{t}\right)$ and $z_{*i}\left(\boldsymbol{f}_{*}^{p}\right)$. For simplicity, we denote $\mu_{t}\left(\boldsymbol{\xi}_{i}\right)=h_{\mu}\left(z_{*i}\left(\boldsymbol{f}_{*}^{t}\right)\right)$ and $p\left(\boldsymbol{\xi}_{i}\right)=h_{p}\left(z_{*i}\left(\boldsymbol{f}_{*}^{p}\right)\right)$, where $h_{\mu}\left(\cdot\right)$ and $h_{p}\left(\cdot\right)$ are some known Lipschitz continuous functions with finite Lipschitz constants. Also, we denote $p_{1}\left(\boldsymbol{\xi}_{i}\right)=p\left(\boldsymbol{\xi}_{i}\right)$ and $p_{0}\left(\boldsymbol{\xi}_{i}\right)=1-p\left(\boldsymbol{\xi}_{i}\right)$.

Given that $\boldsymbol{f}_{*}^{t}$ and $\boldsymbol{f}_{*}^{p}$ may not be uniquely determined in ((ref)), the definitions of $\mu_{t}\left(\boldsymbol{\xi}_{i}\right)$ and $p\left(\boldsymbol{\xi}_{i}\right)$ may also depend on the specific choice of $\boldsymbol{f}_{*}^{t}$ and $\boldsymbol{f}_{*}^{p}$. We omit such dependence in the notations for simplicity, as only the optimal $\boldsymbol{f}_{*}^{t}$'s and $\boldsymbol{f}_{*}^{p}$'s that satisfy Assumption II 6 below are relevant for our discussion and proof.

Then, Theorem (ref), coupled with the conditions on $d_{h*}$, $\beta$, $J$ and $c_{n}$, implies that $\mu_{t}\left(\boldsymbol{\xi}_{i}\right)$ and $p\left(\boldsymbol{\xi}_{i}\right)$ can be estimated with sufficient accuracy for inference. Denote $\hat{\mu}_{t}\left(\boldsymbol{\xi}_{i}\right)=h_{\mu}\left(z_{i}\left(\hat{\boldsymbol{\theta}}_{t}\right)\right)$, $\hat{p}_{1}\left(\boldsymbol{\xi}_{i}\right)=h_{p}\left(z_{i}\left(\hat{\boldsymbol{\theta}}_{p}\right)\right)$, and $\hat{p}_{0}\left(\boldsymbol{\xi}_{i}\right)=1-\hat{p}_{1}\left(\boldsymbol{\xi}_{i}\right)$ as the GNN estimators for $\mu_{t}\left(\boldsymbol{\xi}_{i}\right)$, $p_{1}\left(\boldsymbol{\xi}_{i}\right)$, and $p_{0}\left(\boldsymbol{\xi}_{i}\right)$, respectively, where $\hat{\boldsymbol{\theta}}_{t}$ and $\hat{\boldsymbol{\theta}}_{p}$ are estimated using the following GNN architectures

{

}The estimators $\hat{\boldsymbol{\theta}}_{t}$ and $\hat{\boldsymbol{\theta}}_{p}$ may not be unique, while our discussion holds for any $\hat{\boldsymbol{\theta}}_{t}$ and $\hat{\boldsymbol{\theta}}_{p}$ that solve the minimization problems in ((ref)). Moreover, denote $\boldsymbol{\upsilon}{}_{i}=\left(y_{i},t_{i},\boldsymbol{\xi}_{i}\right)$, and we further introduce the following notations{

}Then the counterfactual policy effect can be identified as $\pi\left(s\right)=\mathbb{E}\left[\zeta_{i}\right]$ under proper conditions, and our GNN-based estimator of $\pi\left(s\right)$ is constructed as{ \[ \hat{\pi}\left(s\right)=\frac{1}{n}\sum_{i=1}^{n}\hat{\zeta}_{i}. \] }In the next subsections, we study the asymptotic distribution of this estimator.

Remark {[}SUTVA assumption{]}: It is important to note that the setup described above relies on the stable unit treatment value assumption (SUTVA), which assumes that an individual's potential outcome only depends on their own treatment assignment. However, a natural extension of this setup allows an individual's potential outcome to also depend on the treatment assignments of other people. This has been extensively studied in the network interference literature (e.g., leung2022causal and savje2021average), where the potential outcome is modeled as $y_{i}(\left\{ t_{j}\right\} _{j\in\left[n\right]})$, instead of $y_{i}\left(t_{i}\right)$.

Although interference is prevalent in many contexts, studying network heterogeneity under SUTVA also holds significant empirical value. In particular, our empirical application studies a substantive policy question about selecting a group of `seed' information recipients, who will then diffuse information throughout social networks via word-of-mouth (Rogers2003diffusion and banerjee2019using). To address this question, it is useful for policymakers to understand the average potential outcome among the chosen seed information recipients under different selection rules.\footnote{This average potential outcome can be learnt once $\pi\left(s\right)$ is known. Specifically, suppose the policymaker considers a selection rule which picks observation $i$ as an information recipient (so $s\left(\boldsymbol{\xi}_{i}\right)=1$) if their $\boldsymbol{\xi}_{i}$ meet certain criteria. Also, for simplicity, assume that $y_{i}\left(0\right)=0$ and the selection probability $Pr\left(s\left(\boldsymbol{\xi}_{i}\right)=1\right)$ is known, as in our empirical application. Then, it follows that $\mathbb{E}\left[y_{i}\left(1\right)\mid s\left(\boldsymbol{\xi}_{i}\right)=1\right]=\frac{\pi\left(s\right)}{Pr\left(s\left(\boldsymbol{\xi}_{i}\right)=1\right)}$, showing that the estimation of $\pi\left(s\right)$ can help the policymaker obtain the estimation of the average potential outcome among the seed information recipients.} This average potential outcome substantially affects the success of the information diffusion process. In the context of BanerjeeChandrasekharDufloJackson2013, for instance, the average potential outcome equates to the probability that a seed household participates in a microfinance program, which largely affects the likelihood that the household subsequently informs its peers about the program. Nevertheless, interference is not necessarily an issue in addressing this policy question. This is because the seed information recipients make their participation decisions before their peers are informed about the program (this is consistent with the first-step model for estimating the characteristic coefficients $\beta$ using the leaders' participation decisions in BanerjeeChandrasekharDufloJackson2013).

Our GNN approach has the potential to contribute to the literature of network interference as well. By including the treatment variable $t_{i}$ as a component of the node attribute vector $\boldsymbol{x}_{i}$, the GNN framework allows an individual's outcome to depend on the treatment assignments of neighbors up to a distance of $L$. Nonetheless, as the empirical focus of our paper is not on network interference, we leave the detailed discussion to future research. $\square$

Assumptions

Assumption II provides the conditions for deriving the asymptotic distribution of the estimator $\hat{\pi}\left(s\right)$ for the average effect of a counterfactual policy $\pi\left(s\right)$.

Assumption II

In particular, the first condition assumes the data is identically distributed, and the second condition requires the individual network heterogeneity is uniformly bounded. Both conditions have been imposed in the derivation of Theorem (ref).

Conditions 3 (a)-(c) assume that the first-stage GNN estimators are well-behaved. This condition can be justified by Theorem (ref) assuming that the prerequisites are fulfilled to ensure the empirical $L_{2}$ convergence rate of the GNN estimators is $o_{p}\left(n^{-1/4}\right)$, as discussed following Theorem (ref). More details for verifying Condition 3 (c) are provided in Appendix (ref).

Conditions 4 and 5 include the unconfoundedness and overlap assumptions, which are standard identification conditions in the treatment effects literature. And the estimated treatment probability is assumed to be bounded inside the interval $\left(0,1\right)$. We assume that the number of layers $L$ in $\boldsymbol{\xi}_{i}$ to fulfill the unconfoundedness condition is known to the researcher, and this $L$ is used throughout this section to construct $\boldsymbol{\xi}_{i}$, to define the network heterogeneities $z_{*i}\left(\boldsymbol{f}_{*}^{t}\right)$ and $z_{*i}\left(\boldsymbol{f}_{*}^{p}\right)$, and to formulate the GNN estimators $z_{i}(\hat{\boldsymbol{\theta}}_{t})$ and $z_{i}(\hat{\boldsymbol{\theta}}_{p})$. We utilize the same $L$ throughout our analysis in this section to simplify our presentation, which however is not a binding constraint.\footnote{We could achieve the same asymptotic distribution in Corollary (ref) by incorporating different numbers of layers across object constructions. In particular, we could allow the counterfactual policy $s\left(\boldsymbol{\xi}_{i,L_{1}}\right)$ to depend on the $L_{1}$-hop neighborhood. And we could define the network heterogeneities $z_{*i}\left(\boldsymbol{f}_{*,L_{2}^{1}}^{1}\right)$, $z_{*i}\left(\boldsymbol{f}_{*,L_{2}^{0}}^{0}\right)$, and $z_{*i}\left(\boldsymbol{f}_{*,L_{3}}^{p}\right)$, as well as the GNN estimators $z_{i}\left(\hat{\boldsymbol{\theta}}_{1,L_{2}^{1}}\right)$, $z_{i}\left(\hat{\boldsymbol{\theta}}_{0,L_{2}^{0}}\right)$, and $z_{i}\left(\hat{\boldsymbol{\theta}}_{p,L_{3}}\right)$, with each depending on the respective $L_{2}^{1}$, $L_{2}^{0}$, and $L_{3}$-hop neighborhoods. Moreover, we could introduce the assumptions that the unconfoundedness condition holds when conditioning on at least the $L_{4}^{t}$-hop neighborhood (i.e., $\mathbb{E}\left[y_{i}\left(t\right)\mid\boldsymbol{\xi}_{i,L},t_{i}\right]=\mathbb{E}\left[y_{i}\left(t\right)\mid\boldsymbol{\xi}_{i,L}\right]$ for $L\geq L_{4}^{t}$ and $t\in\left\{ 0,1\right\} $); the expected potential outcomes depend on at most the $L_{2*}^{t}$-hop neighborhood (i.e., $\mathbb{E}\left[y_{i}\left(t\right)\mid\boldsymbol{\xi}_{i,L}\right]=\mathbb{E}\left[y_{i}\left(t\right)\mid\boldsymbol{\xi}_{i,L_{2*}^{t}}\right]$ for $L\geq L_{2*}^{t}$ and $t\in\left\{ 0,1\right\} $); and the treatment probability depends on at most the $L_{3*}$-hop neighborhood (i.e., $\Pr\left(t_{i}=1\mid\boldsymbol{\xi}_{i,L}\right)=\Pr\left(t_{i}=1\mid\boldsymbol{\xi}_{i,L_{3*}}\right)$ for $L\geq L_{3*}$). Then, the result in Corollary (ref) can be established if the numbers of layers in the GNN estimations satisfy $L_{2}^{t}\geq L_{2*}^{t}$ for $t\in\left\{ 0,1\right\} $ and $L_{3}\geq L_{3*}$. In practice, $L_{2*}^{t}$ and $L_{3*}$ may be unknown, while the authors are currently working on a follow-up paper to formally study the procedure for selecting $L_{2}^{t}$ and $L_{3}$.}

Condition 6 posits that the penultimate layer latent embeddings effectively capture the relevant local neighborhood information in $\boldsymbol{\xi}_{i}$ required to determine the conditional expectations of potential outcomes and the propensity score. Condition 7 assumes that the potential outcome residuals are uncorrelated across observations, conditional on the local neighborhoods and treatment decisions. Condition 8 is a mild regularity condition, assuming that the potential outcomes have finite second moments conditional on the local neighborhoods and treatment decisions.

Condition 9 assumes that $\Sigma_{n}$ is bounded away from zero in the limit. If the data were i.i.d., this condition holds, as $\Sigma_{n}=\text{Var}\left(\zeta_{i}\right)>0$. And $\Sigma_{n}$ would be larger when there is a positive correlation among observations, which is commonly observed in network settings where individuals share common friends. Hence, this condition is easily satisfied.

Condition 10 includes the primitive conditions needed to obtain the central limit theorem for dependent data, which allows $\zeta_{i}$ to be unbounded. If we additionally assume $\zeta_{i}$ is uniformly bounded (i.e., $\max_{i}\left|\zeta_{i}\right|<c_{\zeta}<\infty$ for some fixed constant $c_{\zeta}$), which naturally holds for binary outcome variables, Condition 10 can be replaced with a simpler condition that $\omega_{n}$ grows slower than $\sqrt{n}$ such that $\frac{\left(\omega_{n}\right)^{1-\frac{1}{m}}}{\left(n\right)^{\frac{1}{2}-\frac{1}{m}}}\to0$ for some integer $m$, as discussed in Theorem 2 of janson1988normal. Essentially, Condition 10 requires that the tails of the distribution of $\zeta_{i}$ do not decay too slowly and that $\omega_{n}$ does not grow too fast. In the context of clustered data which we will discuss next, Condition 10 holds true under Assumptions II 2 and 4, along with Assumptions III 2, 3 and 4.

Under Assumption II, we can derive the asymptotic distribution of the estimator $\hat{\pi}\left(s\right)$, which depends on the unknown parameter $\Sigma_{n}$, as presented in Corollary (ref). Then for feasible inference, we introduce some additional conditions, which are suitable for many empirical settings using clustered data, including our empirical context where the data is separated by villages. These conditions are listed under Assumption III.

Assumption III

Condition 1 assumes a clustered structure on the data. We treat the clusters as fixed for each sample size, otherwise the analysis would be overly complicated. So the number of clusters $\bar{c}$ and cluster sizes $\left\{ n_{c}\right\} _{c\in\left[\bar{c}\right]}$ depend on $n$ deterministically. Condition 1 requires that the number of clusters grows with the sample size. This condition suits many empirical datasets with inherent grouping structures, such as students in different schools, employees in various companies, or patients in different hospitals, among others. In our empirical setting, where households are separated by villages, it is also reasonable to allow for free intra-village dependencies while assuming away the possibility of inter-village dependencies. For inference, we apply the same partitions that are used to cluster $\left\{ \zeta_{i}:i=1,\ldots,n\right\} $ to partition their estimators $\left\{ \hat{\zeta}_{i}:i=1,\ldots,n\right\} $, which yields $\left\{ \hat{\zeta}_{cj}:j=1,...,n_{c}\right\} _{c\in\left[\bar{c}\right]}$ used in Corollary (ref).

Condition 2 assumes that each cluster size grows uniformly no faster than the average cluster size. This condition allows the cluster sizes to stay finite. Condition 3 is another primitive condition, which is slightly stronger than assuming that the potential outcomes uniformly have finite second moments. Condition 4 is equivalent to the assumption that $\underset{n\rightarrow\infty}{\lim\inf}\ensuremath{\frac{1}{\bar{c}}}\sum_{c\in\left[\bar{c}\right]}\ensuremath{\left(\frac{n_{c}}{n/\bar{c}}\right)^{2}\mathbb{E}\left[\left(\bar{\zeta}_{c}-\mu\right)^{2}\right]>c_{\sigma}>0}$, where $\bar{\zeta}_{c}\coloneqq\frac{1}{n_{c}}\sum_{j\in\left[n_{c}\right]}\zeta_{cj}$ and $\mu\coloneqq\mathbb{E}\left[\zeta_{i}\right]$, which essentially requires that there is no weak dependence among observations within each cluster.

Asymptotic distribution

In Corollary (ref), we present the asymptotic distribution of the counterfactual policy effect estimator $\hat{\pi}\left(s\right)$ and a feasible estimator of the asymptotic variance for valid inference. The proof of the corollary can be found in Appendix (ref).

Based on the result above, we could construct confidence intervals and conduct hypothesis testing as usual.

Monte Carlo simulations

In this section, we examine the finite sample performance of the average treatment effect ($\tau$) estimator using simulations. The exercise on $\pi\left(s\right)$, the average effect of a counterfactual policy, is very similar and hence omitted in the presentation. In particular, we generate network data, and simulate treatment and outcome variables on the network data. Then, we estimate the GNN models, construct the ATE estimator, and evaluate the coverage of the 95% confidence interval of the ATE estimator. The simulation setup is designed to reflect the conditions of our empirical application.

Specifically, we simulate $\bar{v}=50$ separated networks that we refer to as `villages', with each village containing $n_{v}=200$ observations. For any given node $i$ in a village, the probability of forming a directed edge from another node $j$ within the same village to node $i$ is denoted as $p_{e}$ (i.e., $\text{Pr}\left(d_{ij}=1\right)=p_{e}$). This probability is set to maintain an average of $5$ adjacent neighbors per node, computed as $p_{e}\times(n_{v}-1)=5$. After the initial edge formation, we impose a limit, $c_{n}=10$, on the maximum number of adjacent neighbors a node can have. If a node has more than $c_{n}$ incoming edges (i.e., $\left|\mathcal{N}\left(i\right)\right|\geq c_{n}$), only $c_{n}$ of them will be retained at random.\footnote{We have examined other settings, including those where each node has on average $10$ adjacent neighbors but no more than $40$, which provide similar results and hence are omitted in the presentation.} No edge can be formed across villages. In addition, the covariates, $\boldsymbol{x}_{i}$, are generated i.i.d. uniformly from the $d$-dimensional cube $\left[-1,1\right]^{d}$, where $d$, the number of covariates, is set to $4$. As before, we use $\boldsymbol{\xi}_{i}$ to denote the $L$-hop local neighborhood around node $i$.

We then construct a treatment model and an outcome model. In our empirical application, the treatment decision pertains to whether a household is informed by the experimenter about the microfinance loan, and the outcome decision corresponds to whether a household chooses to participate in the microfinance program if being informed. To emulate the conditions of our empirical application, both the treatment and outcome decisions are binary variables in our simulations.

For the outcome variable, we consider the setting $L=2$, indicating that an observation's decision depends on the covariates of their neighbors and neighbors' neighbors. In particular, we have{

}where we set $d_{h*}=d=4$, $\overline{\boldsymbol{x}}_{*\mathcal{N}\left(i\right)}=\frac{1}{\left|\mathcal{N}\left(i\right)\right|}\sum_{j\in\mathcal{N}\left(i\right)}\boldsymbol{x}_{j}$, and $\overline{\boldsymbol{h}}_{*\mathcal{N}\left(i\right)}=\frac{1}{\left|\mathcal{N}\left(i\right)\right|}\sum_{j\in\mathcal{N}\left(i\right)}\boldsymbol{h}_{*j}$ unless $\mathcal{N}\left(i\right)=\textrm{�}$, otherwise we set $\overline{\boldsymbol{x}}_{*\mathcal{N}\left(i\right)}=\overline{\boldsymbol{h}}_{*\mathcal{N}\left(i\right)}=\mathbf{0}$. We use quadratic models to construct $\boldsymbol{f}_{*}^{(1)}$ and $f_{*}^{(2)}$. Hence, for $l\in\left\{ 1,2\right\} $, we express{

}with $\tilde{\boldsymbol{x}}\in\mathbb{R}^{d^{2}}$ and $\tilde{\boldsymbol{y}}\in\mathbb{R}^{d^{2}}$being the second order interactions of $\boldsymbol{x}$ and $\boldsymbol{y}$, respectively. The outcome model depends on the parameters for the linear terms $\boldsymbol{A}^{(1)},\boldsymbol{B}^{(1)}\in\mathbb{R}^{4\times4}$ and $\boldsymbol{A}^{(2)},\boldsymbol{B}^{(2)}\in\mathbb{R}^{1\times4}$, for which we randomly draw each element from the uniform distribution $\mathcal{U}\left[0.3,0.7\right]$. Also, we randomly draw each element from $\mathcal{U}\left[-0.1,0.1\right]$ for the quadratic-term parameters $\boldsymbol{C}^{(1)},\boldsymbol{D}^{(1)}\in\mathbb{R}^{4\times16}$ and $\boldsymbol{C}^{(2)},\boldsymbol{D}^{(2)}\in\mathbb{R}^{1\times16}$, and from $\mathcal{U}\left[0.3,0.7\right]$ for the constant-term parameters $\boldsymbol{c}^{(1)}\in\mathbb{R}^{4}$, and we set $c^{(2)}=-0.1$. Then, the outcome decision if being treated is generated as $y_{i}\left(1\right)=1\left[z_{*i}\left(\boldsymbol{f}_{*}\right)>\epsilon_{1,i}\right]$, where $\epsilon_{1,i}\sim Logistic\left(0,1\right)$. The outcome decision if not being treated is set as $y_{i}\left(0\right)=0$, which aligns with our empirical application where a household cannot take the microfinance loan if not being informed.

For the treatment model, we consider two scenarios. The first scenario is about random treatment with a constant treatment probability, where the propensity score is set to be $\Pr\left(t_{i}=1\mid\boldsymbol{\xi}_{i}\right)=0.5$, and the treatment decision is generated as $t_{i}=1\left[\mathcal{U}\left[0,1\right]>0.5\right]$. In the second scenario, the treatment decision depends on neighbors up to distance $L=2$, using the same models as for the outcome decision (as shown in ((ref))-((ref))), except that the parameters are different. Specifically, we randomly draw each element of the parameters for linear components, $\boldsymbol{A}^{(l)},\boldsymbol{B}^{(l)}$ for $l\in\left\{ 1,2\right\} $, from $\mathcal{U}\left[0.1,0.2\right]$, for quadratic components, $\boldsymbol{C}^{(l)},\boldsymbol{D}^{(l)}$ for $l\in\left\{ 1,2\right\} $, from $\mathcal{U}\left[-0.1,0.1\right]$, and for constant terms $\boldsymbol{c}^{(1)}$ from $\mathcal{U}\left[0.1,0.2\right]$, and we set $c^{(2)}=-0.1$. Then, the treatment decision is determined as $t_{i}=1\left[z_{*i}\left(\boldsymbol{f}_{*}^{p}\right)>\epsilon_{2,i}\right],$ where $\epsilon_{2,i}\sim Logistic\left(0,1\right)$. After generating the parameters for the treatment and outcome models, these parameters are held fixed across all simulation replications, while the networks, covariates, outcomes, and treatment decisions are regenerated for each simulation replication.

In the simulations, our parameter of interest is the average treatment effect \\ $\tau=\mathbb{E}\left[y_{i}\left(1\right)-y_{i}\left(0\right)\right]=\mathbb{E}\left[y_{i}\left(1\right)\right]$. With the simulated data, we construct the ATE estimator as follows. For the outcome model, we estimate the model with a two-layer GNN using treated individuals only, which provides $\hat{\mu}_{1}\left(\boldsymbol{\xi}_{i}\right)$. And for the treatment model, we use the sample mean of $t_{i}$ for the random scenario, and a two-layer GNN with all individuals for the non-random scenario, which gives $\hat{p}_{1}\left(\boldsymbol{\xi}_{i}\right)$. Then, we compute the ATE estimator as $\hat{\tau}_{n}=\frac{1}{n}\sum_{i\in\left[n\right]}\hat{\zeta}_{i}$, where $\hat{\zeta}_{i}=\frac{\mathbb{I}\left\{ t_{i}=1\right\} }{\hat{p}_{1}\left(\boldsymbol{\xi}_{i}\right)}\left(y_{i}-\hat{\mu}_{1}\left(\boldsymbol{\xi}_{i}\right)\right)+\hat{\mu}_{1}\left(\boldsymbol{\xi}_{i}\right)$, and the estimator of its variance as $\hat{\Sigma}_{n}=\frac{1}{\bar{v}}\ensuremath{\sum_{v\in\left[\bar{v}\right]}}n_{v}(\bar{\hat{\zeta}}_{v}-\bar{\hat{\zeta}})^{2}$, with $\bar{\hat{\zeta}}_{v}$ being the sample mean of $\hat{\zeta}_{i}$ in each village $v$ and $\bar{\hat{\zeta}}$ being the sample mean of $\hat{\zeta}_{i}$ over the full sample.

In terms of the GNN estimation for the outcome and treatment models, the number of hidden neurons for each layer is set to be 8, 16, or 32, which results in nine different architectures for a two-layer GNN estimator. We set the learning rate to be 0.001 and batch size 5. To keep the presentation succinct, we only present results using the same GNN architecture for both the outcome and treatment models. However, we have also experimented with various combinations of architectures, different learning rates (such as 0.005 and 0.01), and different batch sizes (such as 10). These additional results yield very similar findings and hence are omitted from the presentation.

Based on Corollary (ref), we could construct the 95% confidence interval for each simulated sample (i.e., $\hat{\tau}\pm\frac{1.96}{\sqrt{n}}(\hat{\Sigma}_{n})^{1/2}$) and examine if $\tau$ falls within it. We repeat this process for 1000 replications, and report the biases and coverage probabilities of the estimator in Table (ref). Our findings indicate that the ATE estimator exhibits very small biases, and the coverage probabilities are close to the nominal level, for both random and non-random treatment scenarios and across all GNN architectures. The simulation results suggest that the asymptotic distribution in Corollary (ref) is a reasonable approximation to the finite sample dispersion of the estimator.

{}

Empirical application

The estimation of heterogeneous treatment effects is of essential policy relevance. In many empirical contexts, an individual's response to treatment depends not only on their own characteristics but also on their social ties and the characteristics of their neighbors. In this section, we study the empirical context of seeking external financing, which is an important life decision that significantly affects an individual's well-being. In particular, we investigate whether people choose to borrow through a microfinance program when they are informed about this opportunity. From the perspective of policymakers, understanding the probabilities of individual borrowing based on neighborhood information can assist in determining to whom the relevant information should be sent when network data is available.

We utilize the dataset from the paper titled \textquotedbl The Diffusion of Microfinance\textquotedbl by BanerjeeChandrasekharDufloJackson2013. The researchers in this study observe households in 49 isolated Indian villages, with an average village size of 216.69 households. The dataset includes information on household demographics, social networks, and household microfinance participation decisions. In their study, information on microfinance is initially shared with a subset of households in each village, referred to as `leaders'. In total, there are 1262 leaders, with an average of 25.75 leaders per village. The researchers then observe the borrowing decisions made by these leaders after receiving the information. We utilize this dataset to study the average treatment effect of providing loan information, as well as the average outcome of counterfactual policies.

In the realm of microfinance, leaders' borrowing decisions can be influenced not only by their own economic conditions but also by the availability of alternative borrowing options, such as borrowing from their neighbors. Meanwhile, the lending decisions of these neighbors may hinge not only on their own economic conditions, but also on the conditions of others they can in turn borrow from (i.e., the conditions of leaders' neighbors' neighbors). Hence, when modeling leaders' borrowing decisions, we assume that these decisions depend on social ties extending two hops away, and we adopt a two-layer GNN model for estimation. In this estimation, interference is not a primary concern, as leaders are the first to be informed about the loan when others do not have access to such information (this aligns with the model in BanerjeeChandrasekharDufloJackson2013, as we discussed earlier).

For the \textsl{treatment} decisions of the experimenter on disseminating microfinance information, we consider two possibilities: a random assignment with a constant treatment probability and a decision that depends on neighbors up to two hops away, where for the latter we again use a two-layer GNN model for estimation. The authors are working on a subsequent paper that aims to formally determine the relevant neighborhood size for GNN estimations.

In our empirical estimation, we employ the `union relation' from the dataset to represent network edge, which encompasses any social interaction between two households, including activities like visiting each other's homes, going to a temple together, among others. Additionally, we follow the study by BanerjeeChandrasekharDufloJackson2013 and use the following four variables as covariates: number of rooms, number of beds, whether electricity is privately or publicly supplied, and the presence or absence of a latrine. We also include closeness centrality and betweenness centrality as covariates. Our estimation results are presented in the following subsection.

Inference of treatment effects

In Table (ref), we present the estimation of the average treatment effect along with the 95% confidence intervals. As can be seen from the table, the average borrowing probability if being informed about the microfinance program is approximately 24% to 25%, with a relatively tight confidence interval. The ATE estimates are similar across GNN architectures, as well as between the random treatment estimation and the GNN estimation for the treatment decision.

Beyond the ATEs, policymakers may also be interested in understanding the average outcomes of disseminating information to villagers with different network statistics. For example, they may ask: what would the average participation rate be if information is sent to those whose network centrality measure ranks in the top 25th percentile, versus those in the bottom 25th percentile? Specifically, by including household $i$'s centrality measure in $\boldsymbol{\xi}_{i}$, policymakers may compare $\pi\left(s\right)$ where $s\left(\boldsymbol{\xi}_{i}\right)=1$ if $i$'s centrality is above the 75th percentile, with $\pi\left(s\right)$ where $s\left(\boldsymbol{\xi}_{i}\right)=1$ if $i$'s centrality is below the 25th percentile. We focus on two specific centrality measures: closeness centrality and betweenness centrality. Closeness centrality measures the reciprocal of the average of the shortest path distances from one node to all other nodes within the same village, effectively reflecting how quickly a node can reach others. Betweenness centrality, on the other hand, quantifies the number of shortest paths between all pairs of nodes within the same village that pass through a specific node, indicating its role as a `bridge' in the network. Both of these measures are valuable in assessing a node's importance in terms of information dissemination within network structures.

As shown in Table (ref), households with higher centrality tend to exhibit a lower borrowing probability compared to those with lower centrality. The average borrowing probability is approximately $28.0\%$ ($0.070\times4$) to $31.6\%$ ($0.079\times4$) among households within the lowest 25th percentile of closeness centrality, compared to $19.6\%$ ($0.049\times4$) to $20.0\%$ ($0.050\times4$) among households within the top 25th percentile of the same measure. Similar patterns are found for betweenness centrality. We have also experimented with degree centrality and eigenvector centrality, both of which produced similar findings and are hence omitted for brevity. A plausible explanation for these findings could be that households with higher centrality tend to have more avenues for borrowing within their social networks, thereby reducing their need to participate in the microfinance program.

According to the study by BanerjeeChandrasekharDufloJackson2013, a selected group of households (i.e., leaders) is informed about the microfinance program and subsequently decides whether to participate. These leaders then have a certain probability of informing their acquaintances about the program, thereby facilitating the spread of knowledge about the microfinance program. Importantly, the likelihood of leaders disseminating this information depends on their own participation decisions. As demonstrated empirically in Table 1 of the paper by BanerjeeChandrasekharDufloJackson2013, the probability of a participating leader disseminating the information is 7-10 times higher than that of a non-participating leader.

Given budget constraints, policymakers often face the challenge of deciding who should initially receive the program information, as it would be prohibitively costly to inform everyone. The aim is to distribute this information as quickly and as widely as possible through social networks to enhance overall program participation. However, solely targeting households based on their centrality measure might not be the most effective approach, as the participation probability of the targeted households also plays a significant role in information diffusion. This consideration prompts our exploration of more balanced targeting strategies in the following subsection, aiming to identify and prioritize individuals who exhibit not only high centrality but also a strong willingness to participate.

Individual targeting on networks

As discussed previously, the success of the program hinges on a judicious selection of leaders. The process must balance the leaders' likelihood of participating with the strength of their social ties. The program is likely to flourish if the chosen leaders not only are inclined to participate but also maintain strong connections with uninformed households in their social networks. This is an essential crux of network diffusion problems. To address this, we propose an improved approach to leader selection that outperforms the current method.

We use the GNN estimation to predict the participation probabilities of households based on their local neighborhoods if they were informed about the program. We denote the predicted participation probability for household $i$ as $\hat{p}_{i}$. To measure the strength of social ties for each household $i$ in the network, denoted as $c_{i}$, we use either closeness centrality or betweenness centrality separately. We have also explored degree centrality and eigenvector centrality, which yielded similar results and hence are omitted in the presentation.

To select a group of households that balance the participation probability $\hat{p}_{i}$ and centrality $c_{i}$, we use a relative weight $\omega\in\left[0,1\right]$ to capture the balance between these two aspects, and we compute a weighted score for each household as follows{ \[ s_{i}=\left(1-\omega\right)\hat{p}_{i}+\omega c_{i}. \] }As $\omega$ approaches 0, the score places more emphasis on participation probability. Conversely, when $\omega$ approaches 1, the score assigns greater emphasis to centrality. Fixing a relative weight $\omega$, we can select the top $k$ households based on the weighted score. This gives us a group of the $k$ best leaders that balance the leader participation and centrality at this particular weight $\omega$. For our exercises, we set $k$ to equal 1262, which is the actual number of leaders selected by the researchers in BanerjeeChandrasekharDufloJackson2013.

By varying $\omega$, we can find different groups of $k$ best leaders and compute their average participation probabilities and centralities. This enables us to trace out the Pareto frontier that the program can achieve between leader participation rate and centrality, as indicated by the blue curves in Figure (ref). The Pareto frontier illustrates potential improvements compared to the selected leaders under the current strategy, as marked by the orange star in Figure (ref). Regarding centrality and participation probability, it is possible to realize improvements equivalent to 1.5 standard deviations in betweenness centrality, 0.94 standard deviations in closeness centrality, or 1.83 standard deviations in leader participation probability alone. Also, there is a large room for joint Pareto improvements in both participation probability and network centrality, relative to the current leader selection rule. Our findings present valuable insights on the enhancement of information diffusion through social networks.

Conclusion

In this paper, we present a novel application of graph neural networks for modeling and estimating network heterogeneity, focusing on the empirical context where individual outcomes or decisions depend on local neighborhood surroundings. We provide theoretical justifications for the GNN estimator and demonstrate its usage in causal inference with heterogeneous treatment effects. Empirically, in the context of a microfinance program, we show that the estimator can be applied to evaluate various treatment effects and to enhance information targeting on networks. We believe the GNN estimator holds substantial value in economics and social sciences, leveraging the potency of versatile neural network architectures and abundant network data. Adapting the GNN model to accommodate diverse empirical needs in different economic settings opens up an exciting avenue for future research.