Difference-in-Differences using Double Negative Controls and Graph Neural Networks for Unmeasured Network Confounding

Difference-in-Differences using Double Negative Controls and Graph Neural Networks for Unmeasured Network Confounding

\def\spacingset#1{ {#1}} \spacingset{1}

\if11 \fi

\if01 {

} \fi

{\it Keywords:} causal inference; network interference; double robustness; high-dimensional data; approximate neighborhood interference

\spacingset{1.8}

Introduction

Estimating causal effects under network interference has become increasingly important in various fields. For example, scholars have explored the diffusion of microcredit within households he2024measuring, as well as the effect of anti-conflict interventions on adolescent social norms paluck2016changing. In the field of public health, researchers have not only examined the transmission effects of infectious diseases halloran1995causal,morozova2018risk, but also investigated how health behaviors “contagion” among populations christakis2013social. When individuals form connections through a network, identification and estimation of causal effects have been challenging for two reasons. The first issue is that the conventional potential outcome framework rubin1980randomization, which relies on the Stable Unit Treatment Value Assumption (SUTVA), is inappropriate. A large body of literature has focused on randomized controlled trials toulis2013estimation,athey2018exact,li2022random. A common method to relax the SUTVA is to suppose a low-dimensional function, which is called the effective treatment manski2013identification or exposure mapping 10.1214/16-AOAS1005. This function serves as a sufficient statistic for spillover effects, meaning that others' interventions influence an individual's outcome only through it. The exposure mapping is useful for summarizing potentially complex spillover effects, yet there is an inherent challenge in determining the “right” functional form. Consequently, several recent studies have explored the conditions under which meaningful causal parameters can be estimated even when interference is unknown savje2021average,leung2022,savje2024causal,hoshino2024causal. Beyond randomized controlled trials, a growing number of studies have considered observational causal inference in settings with network interference van2014causal,forastiere2021identification,tchetgen2021auto,leunggnn,leung2023network,ogburn2024causal. These approaches typically rely on the unconfoundedness assumption, often requiring conditioning on the covariates of an individual and their neighbors Liu2019, or on the covariates of all individuals in the network leunggnn. However, a key limitation of these methods is their reliance on the assumption of no unmeasured network confounding. This assumption is frequently violated in observational settings, complicating the identification and estimation of causal effects. Such unmeasured confounding generally arises from two sources: the homophily bias (connections driven by unobserved characteristics) and the contextual confounding (peers sharing unobserved contextual factors) manski1993identification,vanderweele2013social. The two aforementioned issues can be summarized as the problem of unmeasured network confounding in observational studies. In this paper, we employ the double negative controls (DNC) and graph neural networks (GNNs) to address this challenge. The DNC specifies a negative control outcome---an outcome variable known not to be causally affected by the treatment---and a negative control exposure, which is a treatment variable that does not causally affect the outcome lipsitch2010negative. In recent years, a growing body of research has used DNC as proxy variables to identify causal effects in the presence of unmeasured confounding, an approach known as proximal causal inference miao2018identifying,shi2020multiply,cui. These methods all rely on the SUTVA and do not involve network data. Concurrently, we introduce GNNs to address the high-dimensional estimation challenges posed by network confounding in our DNC framework. leunggnn has demonstrated that GNNs are particularly well-suited for network-structured data, capable of capturing one's own high-dimensional covariates, neighbors' high-dimensional covariates, and their complex network relationships, thereby enabling effective estimation. xu proposed a DID framework for network interference, but the identification strategy relies on the selection on observables assumption and thus cannot address the issue of unmeasured network confounding. Our contribution is to propose a general Difference-in-Differences (DID) framework that integrates DNC and GNNs. Based on our framework, researchers are able to identify and estimate the direct and indirect causal effects of the treatment in settings with unmeasured network confounding. To achieve this, our identification strategy addresses the limitations of the standard parallel trends assumption by instead positing a latent parallel trends (cf. Assumption 2.1). Furthermore, this assumption is weaker than the network version of unconfoundedness leunggnn,egami and forms the basis for our DNC approach. The nonparametric identification of the effects of interest is derived by leveraging DNC that are associated with unmeasured confounders. Specifically, we incorporate DNC via either an outcome confounding bridge function or a treatment confounding bridge function, which are the network-adapted versions of those studied in cui. Identification is achieved as long as at least one of these two bridge functions satisfies the identification assumptions. Moreover, the semiparametric proximal causal inference framework permits nonparametric estimation of these quantities. Building on this identification strategy, we propose the doubly robust DID estimators. drdid proposed doubly robust DID estimators for independent data, whereas our method accommodates settings with network confounding. Combining the generalized method of moments (GMM) hansen1982large and GNNs, we develop an approach to handle the high-dimensional nonparametric terms within the bridge functions. Subsequently, the asymptotic normality of this estimator is established under assumptions of $\psi$-network dependence kojevnikov2021limit and approximate neighborhood interference (ANI) leung2022. Finally, simulations demonstrate the finite-sample performance of the estimator, and we apply the method to evaluate the impact of the green credit policy on green innovation. The rest of the article is organized as follows. Section (ref) details our setup, identification challenges, and our identification strategy using DNC. Section (ref) presents the GNN-based estimation and inference. Section (ref) evaluates our method via simulations. Section (ref) provides an empirical application, and Section (ref) concludes.

Setup and Identification

Setup

Consider a set of units $N_{n}$ = $\{ 1, 2, \dots, n \}$. The units form an undirected network represented by the $n$ $\times$ $n$ symmetric adjacency matrix $\mathbf{A}$ = $(A_{ij})_{i,j\in N_{n}}$, where $A_{ij}$ $\in$ $\{ 0, 1 \}$ indicates whether or not $i$ and $j$ are connected. We assume that there are no self-links so that $A_{ii}$ = 0 for all $i$ $\in$ $N_{n}$. Define $\mathcal{A}_n = \{0,1\}^{n \times n}$ as the space of all binary adjacency matrices. Let $Y_{it}$ $\in$ $\mathbb{R}$ and $D_{it}$ $\in$ $\{ 0, 1 \}$ denote the observed outcome and treatment status, respectively, for individual $i$ at time $t$, where $t$ $\in$ $\{0,1\}$. Because individuals are exposed to treatment only at $t$ = 1, we have $D_{i0}$ = 0 for all $i$. To reduce notation, we define $D_{i}$ $\equiv$ $D_{i1}$. Denote the $n$-dimensional vector of realized treatments as $\mathbf{D}$ = $(D_{i})_{i \in N_{n}}$, with the support $\mathcal{D}_{n}$ = $\{0,1\}^{n}$. For each $\mathbf{d}$ $\in$ $\mathcal{D}_{n}$, let $Y_{it}(\mathbf{d})$ denote unit $i$'s potential outcome at time $t$ under treatment assignment $\mathbf{D}$ = $\mathbf{d}$. By construction, we have $Y_{it}$ = $Y_{it}(\mathbf{D})$. Denoting $\mathbf{d_{-i}}$ = $(d_{k})_{k\ne i}$, we write the potential outcome of unit $i$ at time $t$ as $Y_{it}(d_{i}, \mathbf{d_{-i}})$, given $D_{i}$ = $d_{i}$ and $\mathbf{D_{-i}}$ = $\mathbf{d_{-i}}$. Because only one realization from $(Y_{it}(\mathbf{d}))_{\mathbf{d} \in \mathcal{D}_{n}}$ is observable for each unit, it is generally impossible to define identifiable causal estimands without introducing restrictions hoshino2024causal. Accordingly, we define a pre-specified function $G$: $N_{n}$ $\times$ $\{ 0, 1 \}^{n-1}$ $\times$ $\mathcal{A}_{n}$ $\to$ $\mathcal{G}$, where $\mathcal{G}$ $\subset$ $\mathbb{R}^{dim(G)}$ is a set that does not depend on $i$ and $n$, and $dim(G)$ is a fixed positive integer. The exposure mapping $G$ defines the exposure realization for unit $i$ as $G_{i} = G(i, \mathbf{D_{-i}}, \mathbf{A})$, which maps the treatment statuses of other units ($\mathbf{D_{-i}}$) and the network structure ($\mathbf{A}$) to a specific exposure value. Following savje2024causal, we employ $G$ to define causal effects without necessarily assuming it to fully capture the complete causal structure. For example, $G_{i}$ = $\sum_{j=1}^{n}A_{ij}D_{j}$. If $G$ is correctly specified, it implies that individual i's outcome depends only on its own treatment status and the number of treated neighbors. However, even under misspecification, $G$ remains a useful construct for defining causal estimands of interest, such as spillover effects arising from variation in the number of treated neighbors. We focus on two types of estimands defined by exposure mappings. The first estimand is the average direct effect on the treated (ADT) at exposure level $g \in \mathcal{G}$, defined as:

For instance, if $G_{i}$ is defined as the number of treated neighbors ($G_{i} = \sum_{j=1}^{n}A_{ij}D_{j}$), then the $\tau_{ADT}$(3) represents the conditional average difference in potential outcomes between being treated and untreated for individuals who are treated and have exactly three treated neighbors. We now examine the average indirect effect on the treated (AIT) estimands. Let $\ell_{\mathbf{A}}(i,j)$ denote the path distance between units $i$ and $j$, defined as the length of the shortest path connecting them. By convention, we set $\ell_{\mathbf{A}}(i,j)$ = $\infty$ when no path exists between $i$ and $j$ in $\mathbf{A}$ and 0 if $i$ = $j$. Suppose that each $G_{i}$ depends only on unit $j$'s such that 1 $\le$ $\ell_{\mathbf{A}}(i,j)$ $\le$ $K$ with some constant $K$ $\ge$ 1 (see Assumption 4.4). Based on this, we define the interference graph $\mathbf{E}$ = $(E_{ij})_{i,j\in N_{n}}$, where $E_{ij}$ = $\mathbf{1}\{ 1 \le \ell_{\mathbf{A}}(i,j) \le K \}$. For each $i$ $\in$ $S_{n}$, let interference set $\mathcal{E} _{i}$ = $\{ j \in N_{n}: E_{ij} = 1 \}$. The AIT estimand is given by

For example, when $K$ = 2, the $\tau_{AIT}$ measures the average indirect causal effect of a treated unit's intervention status on its direct neighbors (distance 1) and neighbors of neighbors (distance 2) through network paths. Units' potential outcome may not share a common expectation, so the causal estimand is explicitly written as the network average of individual level causal effects. Due to space limitations, our discussion primarily focuses on ADT, while a detailed discussion of AIT is provided in Appendix C.

Identification challenge

The DID is one of the most popular methods in the social sciences for estimating causal effects in observational studies roth2023s. The application of the DID to identify the effect of interest under cross-unit interference presents two main challenges. The first challenge is that when SUTVA is relaxed, the conditional parallel trends assumption does not suffice to identify the causal effect. A natural approach to adapting the conditional parallel trends assumption for identifying the ADT under cross-unit interference can be expressed as:

where $\mathbf{X}$ = $(X_{i})_{i=1}^{n}$ represents the matrix of all units' observed covariates, and $X_{i}$ $\in$ $\mathbb{R}^{d}$. However, due to the presence of unmeasured confounding $\mathbf{U}$ = $(U_{i})_{i=1}^{n}$, this assumption generally fails to hold. This assumption is similar to Assumption 3 in xu. Figure 1 shows a causal-directed acyclic graph (DAG). Without loss of generality, we assume there are only two samples, $\Delta Y_{i} = Y_{i1} - Y_{i0}$ for $i \in \{1, 2\}$. If $D_{1} = 1$ and $D_{2} = 0$, then $D_{1} \to \Delta Y_{1}$ and $D_{1} \to \Delta Y_{2}$ denote the ADT and the AIT, respectively. When equation (3) holds, there exists an unblocked back-door path $D_{1}$ $\gets$ $\mathbf{U}$ $\to$ $\Delta Y_{1}$. Consequently, even in the simple setting of dyadic data, identification of the ADT requires additional assumptions. As the identification challenges for the AIT are analogous, we omit a separate discussion here.

Identification with double negative controls

In this section, we consider an alternative identification strategy that leverages auxiliary variables known as negative controls to address unmeasured confounding. While measuring all confounding factors would ensure the standard conditional parallel trends assumption (3) holds, this is often infeasible. We therefore posit the following latent parallel trends assumption.

This is similar to the latent ignorability of egami. Assumption 2.1 states that $\mathbf{U}$ and $\mathbf{X}$ suffice to account for confounding, whereas $\mathbf{X}$ alone may not. This assumption is often plausible as there is no direct restriction on the nature of the latent characteristic $\mathbf{U}$. Additionally, we impose the no anticipation assumption, stipulating that units do not alter their behavior in period $0$ in expectation of future treatment assignment.

This assumption is widely adopted in the studies on the DID models with interference butts2021difference,deuchert2019direct. However, Assumptions 2.1-2.2 are not sufficient for identification due to the unmeasured variables $\mathbf{U}$. For identification under the unobserved confounding, we incorporate two observed auxiliary variables, NCE $Z$, and NCO $W$, which satisfy the following conditions.

Assumption 2.3(a) defines the NCO property, which states that $W_{i}$ is conditionally independent of the treatment $D_{i}$ given the unmeasured confounders $\mathbf{U}$, observed covariates $\mathbf{X}$, network structure $\mathbf{A}$, and exposure mapping $G_{i}$. Assumption 2.3(b) defines the NCE property, which states that $Z_{i}$ is conditionally independent of both the outcome $\Delta Y_{i}$ and the NCO $W_{i}$, given $\mathbf{U}, \mathbf{X}, \mathbf{A}, G_{i}$, and $D_{i}$. Identifying the ADT using NCO W and NCE Z hinges on a key assumption within the proximal causal inference framework miao2018identifying: the existence of a confounding bridge function, as specified in Assumption 2.4.

Assumption 2.4 states that the confounding effect of $\mathbf{U}$ on outcome $\Delta Y_{i}$ is equal to the confounding effect of $\mathbf{U}$ on outcome confounding bridge function $h_1$($W_{i}$, $\mathbf{X}$, $\mathbf{A}$), a transformation of $W_{i}$. Since the first term in the ADT contrast is point-identified under the consistency assumption alone, we only requires invoking the confounding bridge function for the counterfactual outcome under the control group. Equation (4) is formally a Fredholm integral equation of the first kind kress1989linear. The existence of a solution to this equation, as detailed in Appendix A.1 of egami, demands that $W$ be sufficiently informative for $U$, in addition to a set of regularity conditions.

To achieve identification under unobserved confounding, we next employ the NCE $Z$.

This assumption states that $Z$ contains sufficient information about $W$, which is crucial for identifying the outcome confounding bridge function $h_1$. It is a well-known technical condition in the study of sufficiency in statistical inference, called the completeness condition. In practice, selecting NCO and NCE that satisfy Assumption 2.5 is crucial for effect estimation. egami discuss the impacts on estimation when the assumption is violated. To satisfy the completeness condition for categorical variables, the number of categories in NCE needs to be at least as large as the number of categories in NCO. For continuous variables, the number of NCEs needs to be at least as large as the number of NCOs.

The selection of valid negative controls in practice is crucial. Figure 2 illustrates an example of how to choose negative controls that satisfy Assumption 2.2. We consider a case with four samples where the effect of interest is $D_{1} \to \Delta Y_{1}$. In this DAG framework, the treatment assigned to each unit exerts direct effects only on itself and its immediate neighbors. $D_{2}$ and the two variables $\{D_{3}, D_{4}\}$ satisfy the NCO and NCE conditions, respectively. This negative control selection approach is applicable when there are no direct causal relationships between individual treatments, and requires that the number of first-order neighbors (distance $s$ = 1) is substantially smaller than the number of individuals at distance $s$ $\ge$ 2.

Under the stated assumptions, we establish the nonparametric identification of the ADT.

Notably, we identify the ADT without imposing any parametric restriction on the confounding bridge function, and the proof of Theorem (ref) can be find in Appendix A.1.

Doubly robust difference-in-differences estimands

In this section, we examine an alternative identification approach cui and subsequently derive the doubly robust DID estimands.

Analogous to Assumptions 2.4 and 2.5, Assumptions 2.6 and 2.7 allow us to establish another non-parametric identification of the ADT.

Theorem (ref) provides a new identification result for the ADT, and its proof can be found in Appendix A.2. To estimate causal effects, a common approach involves specifying a parametric or semi-parametric function for the confounding bridge function $h_1(W_{i}, \mathbf{X}, \mathbf{A}; \gamma_1)$ with $\gamma_1$, or $q_1(Z_{i}, \mathbf{X}, \mathbf{A}; \gamma_2)$ with $\gamma_2$ miao2018identifying,egami,cui. Then, based on the conditional moment restriction in Theorem 1 (for $h_1$) or Theorem 2 (for $q_1$), one can treat the conditioning set as a high-dimensional vector function and estimate $\gamma_1$ or $\gamma_2$ via GMM method. However, this method relies on the correct specification of the outcome confounding bridge function and the treatment confounding bridge function. If either the function $h_1$ or $q_1$ are misspecified, the plug-in estimator may be severely biased. To address this challenge, we propose the doubly robust DID estimand tailored for unmeasured network confounding. Here, the doubly robust property implies that the resulting estimand identifies the the causal effect of interest, as long as at least one of confounding bridge functions is correctly specified. Let $h^*_1(W_{i}, \mathbf{X}, \mathbf{A})$ and $q^*_1(Z_{i}, \mathbf{X}, \mathbf{A})$ be arbitrary models for the true, unknown function $h_1(W_{i}, \mathbf{X}, \mathbf{A})$ and $q_1(Z_{i}, \mathbf{X}, \mathbf{A})$ for $i$ $\in$ $N_n$. In order to describe our proposed doubly robust approach, consider the following two models:

The doubly robust DID estimand can be expressed as:

The following Proposition (ref) shows that our proposed doubly robust estimand recovers the ADT provided that at least one of models $\mathcal{M}_1$ and $\mathcal{M}_2$ is correctly specified.

The proof can be find in Appendix A.3.

Estimation and Inference

In this section, we discuss the estimation and inference for ADT while accounting for network-dependent in observational studies.

Doubly robust difference-in-differences estimators using graph neural networks

Our estimation procedure can be decomposed into two steps. In the first step, GNNs are employed to estimate the working nuisance models $h^*_1(W_{i}, \mathbf{X}, \mathbf{A})$ and $q^*_1(Z_{i}, \mathbf{X}, \mathbf{A})$, $i$ $\in$ $N_n$, $g$ $\in$ $\mathcal{G}$. In the second step, one plugs the fitted values of the working nuisance models into the sample analogue of $\tau^{dr}_{ADT}(g)$. The standard GNNs architecture comprises $L$ layers, each containing $n$ computational units known as neurons. Each neuron, denoted as $h_i^{(l)}$ for the $i$-th neuron in layer $l$, implements a parameterized, vector-valued function. The connectivity between layers is governed by adjacency matrix $\mathbf{A}$ via a “message-passing” mechanism, for layers $l$ = 1, $\dots$, $L$,

where $\Phi_{0l}$, $\Phi_{1l}$ are parameterized, vector-valued functions, and $h_i^{(0)}$ = $X_i$ denotes the initial node, which contains no inherent network information. In subsequent layers, unit $i$'s node embedding is a function of its $1$-neighborhood's embeddings in the previous layer and therefore incorporates increasingly more network information as $l$ increases. For estimating the confounding bridge functions, we utilize a strategy that combines GNNs with the GMM method. Let $\mathcal{F}_{GNN}(L)$ denote the set of all GNNs with $L$ layers ranging over all possible functions $\Phi_{0l}$, $\Phi_{1l}$ for $l$ = 1, $\dots$, $L$ within some function class. For any $f$ $\in$ $\mathcal{F}_{GNN}(L)$, let $f(i, \mathbf{X}, \mathbf{A})$ denote its $i$th component, which corresponds to $h_i^{(L)}$. To allow for flexible non-linear interactions while strictly satisfying the identification condition that bridge functions depend on unit-specific negative controls, we employ a late-fusion neural network architecture. Suppose both functions $h^*_1(W_{i}, \mathbf{X}, \mathbf{A})$ and $q^*_1(Z_{i}, \mathbf{X}, \mathbf{A})$ take the following form:

where $\Psi_h$ and $\Psi_q$ are multilayer perceptrons (MLPs) that take the concatenation of the negative control variable and the GNN-learned embedding $f(i, \mathbf{X}, \mathbf{A})$ as input. This formulation replaces the restrictive linear specification with a general non-parametric form. Based on equation (5), a moment function is defined for each $i$ as:

where $(\mathbf{A}\mathbf{X})_{i\cdot} \in \mathbb{R}^d$ denotes the $i$-th row of the matrix product $\mathbf{A}\mathbf{X}$, which is the vector sum of the feature vectors of all individuals directly connected to individual $i$. Then, the estimator of $h^*_1$ is

where $\bar{m}(h^*_1)=1/n\sum_{i\in N_n}m_h(\Delta Y_i,W_i,D_i,G_i,Z_{i}, \mathbf{X}, \mathbf{A})$, and $\Omega$ is a user-specified positive-definite weight matrix. Next, $q^*_1(Z_{i}, \mathbf{X}, \mathbf{A})$ is estimated based on (6). Similarly, define a moment function for each $i$ as

Then, the estimator of $q^*_1$ is

where $\bar{m}(q^*_1)=1/n\sum_{i\in N_n}m_q(W_i,D_i,G_i,Z_{i}, \mathbf{X}, \mathbf{A})$. To obtain an efficient estimator, we use the two-step GMM to get the optimal estimator hansen1982large. In the first step, we choose an identity matrix as $\Omega$ or some other positive-definite matrix, and compute the preliminary estimate $\hat{h}^*_1(1)$. In the second step, we compute $\hat{\Lambda}$ by

Then, the final estimate is obtained using the weighting matrix $\Omega$ = $\hat{\Lambda}^{-1}$. The $q^*_1$ is estimated in a similar manner.

The number of layers $L$ in a GNNs determines its receptive field, meaning each i's estimate is constructed using information exclusively from its $L$-neighborhood. As discussed in leunggnn, the ANI assumption (see Assumption 3.4) is analogous to the approximate sparsity assumption in lasso literature, demonstrating that selecting a small $L$ and making estimation of the nuisance functions feasible. Finally, the doubly robust estimator of ADT is given by

where $\hat{\mathrm{E}}[\mathbf{1}_i(g,1)]$ = $1/n\sum_{i\in N_n}\mathbf{1}_i(g,1)$, $h^*_1(W_{i}, \mathbf{X}, \mathbf{A})$ and $q^*_1(Z_{i}, \mathbf{X}, \mathbf{A})$ are constructed based on the aforementioned GNNs estimators. The following assumption imposes restrictions on the convergence rate of the GNNs estimators.

These are similar with the standard conditions farrell2015robust for machine learners. leunggnn examines the validity of this assumption for GNNs estimators and chen2024causal provides the rate of convergence for artificial neural networks under certain regularization conditions. Assumption 3.1 (a) is a mild consistency requirement for the GNNs estimators. Assumption 3.1 (b) requires an explicit rate for the product of errors. Thus, if one function is relatively easy to estimate, Assumption 3.1 (b) can still be satisfied even when the other does not converge at the $n^{-1/4}$. Assumption 3.1 (c) is similar with Assumption 3 (c) of farrell2015robust, it constrains the convergence rate of the product between the GNNs estimators' error and some terms whose expectations are zero.

Inference

We next examine the convergence properties of the ADT estimator and impose the following conditions.

Assumption 3.2 is standard condition and can be generalized to uniformly bounded moments. Assumption 3.3 ensures that the set used for estimation is non-empty. For a nonnegative integer $s$ $\ge$ 0, let $N_{\mathbf{A}}(i,s)$ = $\{j \in N_n: \ell_{\mathbf{A}}(i,j) \le s \}$ denote the set of units within a distance $s$ from unit $i$, referred to as the $s$-neighborhood of unit $i$ (Note that $i$ $\in$ $N_{\mathbf{A}}(i,s)$). Define $\mathbf{d}_{N_{\mathbf{A}}(i,s)}$ = ($d_j: j \in N_{\mathbf{A}}(i,s)$) and $\mathbf{A}_{N_{\mathbf{A}}(i,s)}$ = ($A_{kl}: k,l \in N_{\mathbf{A}}(i,s)$), respectively the subvector of $\mathbf{d}$ and subnetwork of $\mathbf{A}$ on $N_{\mathbf{A}}(i,s)$. Additionally, let $N^c_{\mathbf{A}}(i,s)$ = $N_n \texttt{\textbackslash} N_{\mathbf{A}}(i,s)$ denote the set of units who are more than distance $s$ away from $i$.

The Assumption 3.4 (a) states that the exposure mapping of each unit depends only on the unit's own K-neighborhood. This is a weak restriction on $G$ satisfied by most exposure mappings of interest in the literature leung2022. And Assumption 3.4(b) implies that $D_i$ is unrelated to the treatment of individuals beyond its $K$-neighborhood. Below we introduce the ANI assumption. Let $\mathbf{D'}$ be an independent copy of $\mathbf{D}$, define $\mathbf{D}_i^{(s)}$ = $(\mathbf{D}_{N_{\mathbf{A}}(i,s)}, \mathbf{D}'_{N^c_{\mathbf{A}}(i,s)})$ obtained by concatenating the subvector of $\mathbf{D}$ on $N_{\mathbf{A}}(i,s)$ and that of $\mathbf{D}'$ on $N_n\texttt{\textbackslash}N_{\mathbf{A}}(i,s)$. Finally, let

This requires interference from distant alters to be negligible for large distances. The ANI condition is substantially less restrictive than the conventional clustered interference assumption, which imposes the stringent requirement that $\theta_{n,s}^{ADT}$ =0 for some finite cutoff distance $s$. Given our analysis involves a network of arbitrarily interconnected units, the independently and identically distributed assumption becomes untenable. To derive the asymptotic distribution, we employ the central limit theorem (CLT) for $\Psi$-network dependence established by kojevnikov2021limit.

This work extends temporal $\Psi$-dependence doukhan1999new to network settings by replacing temporal distance with shortest-path distance. The $\Psi$-network dependence requires that the covariance between any two observation sets $\mathbf{C}_H$ and $\mathbf{C}_{H'}$ decays to zero as the distance between them increases. Define

and define $\sigma_{n}^{ADT}$ = Var$\left[n^{-1/2}\sum_{i \in N_n}\left(C_i-\tau^{dr}_{ADT}(g)\right) \right]$. Let $K$ be the constant in Assumption 3.4, $\lfloor s \rfloor$ be $s$ rounded down to the nearest integer, $i^{*}(\cdot)$ be the identity function $x \mapsto x$ on $\mathbb{R}$, and Lip($f$) be the Lipschitz constant and $\|f\|_{\infty}$ = $\sup |f|$.

The proof of Theorem (ref) can be found in Appendix. Theorem (ref) shows that $\{C_i\}_{i=1}^n$ is conditionally $\Psi$-dependent, and we can apply results due to kojevnikov2021limit to show that $\hat{\tau}^{mr}_{ADT}(g)$ is asymptotically normal. Let $N^{\partial}_{\mathbf{A}}(i,s)$ = $\{j \in N_n: \ell_{\mathbf{A}}(i,j) =s\}$ be the subset of $N_n$ that are exactly at distance $s$ from unit $i$ $\in$ $N_n$. Then, define $M^{\partial}_{N_n}(s,k)$ = $n^{-1}\sum_{i\in N_n}|N^{\partial}_{\mathbf{A}}(i,s)|^k$, which measures the denseness of $\mathbf{A}$ restricted on $N_n$. When $k$ = 1, we denote $M^{\partial}_{N_n}(s)$ = $M^{\partial}_{N_n}(s,1)$. Furthermore, let $\Delta_{N_n}(s, m; k)$ = $n^{-1}\sum_{i \in N_n}max_{j \in N^{\partial}_{\mathbf{A}}(i,s)}|N_{\mathbf{A}}(i,m)\backslash N_{\mathbf{A}}(j,s-1)|^k$, where $N_{\mathbf{A}}(j,s-1)$ = $\varnothing$ if $s$ = 0. This represents the $k$th sample moment of the maximum number (across all $j$ at distance s from $i$) of units who are within distance $m$ from $i$ but at least distance $s$ apart from $j$. In addition, we define $c_{N_n}(s,m;k)$ = $inf_{\alpha > 1}$$[\Delta_{N_n}(s, m; k\alpha)]^{\frac{1}{\alpha}}$$[M^{\partial}_{N_n}(s,\alpha/(\alpha-1))]^{1-\frac{1}{\alpha}}$. This quantity measures the denseness of the network, which plays an important role in establishing the CLT hoshino2024causal.

Assumption 3.6(a) rules out that there are a non-negligible proportion of units whose $2K$ neighborhoods may grow to infinity as $n$ increases. If one assumes that each individual can hold only a limited number of interacting partners, Assumption 3.6(a) is satisfied with $M_{S_n}^\partial (s, k) < \infty$ for all $s, k < \infty$. However, for example, it is violated if the network is a complete graph. The first two terms in Assumption 3.6(b) correspond to Assumption 3.4 of kojevnikov2021limit, which they utilize to establish a CLT. The third is similar with leunggnn and used to asymptotically linearize our robust estimator under network dependence.

The proof of Theorem (ref) can be find in Appendix A. The next result characterizes the asymptotic properties of $\hat{\sigma}_n^{ADT}$. We consider inference methods based on network HAC estimation and demonstrate that the HAC estimator exhibits asymptotic biases. This bias arises from the inability to estimate the heterogeneous means that appear in the asymptotic variances in Theorem (ref). This is a well-known issue in the design-based uncertainty framework imbens2015causal. Define

Part (a) provides a mild strengthening of Assumption 3.1, as all nuisance functions are estimated nonparametrically. Since it does not require uniform convergence, it is easier to verify for machine learning estimators. Parts (b)-(e) restrict both the network structure and the rate of divergence of $b_n$ in a similar manner to Assumption 7 of leunggnn and Assumption 4.1 of kojevnikov2021limit. Define

where

Similarly, define

where $\tilde{\tau}_{ADT,i}(g)$ = $C_i -\mathrm{E}[C_i]$, and

The proof of Theorem (ref) can be find in Appendix A.6. In this article, we apply the bandwidth

where $\lceil \cdot \rceil$ rounds up to the nearest integer, $\delta(\mathbf{A})$ = $n^{-1}\sum_{i,j} A_{ij}$ is the average degree, and $\mathcal{L}(\mathbf{A})$ = $(n(n-1))^{-1}\sum_{i\ne j} \ell_{\mathbf{A}}(i,j)$ is the average path length, which means the average over all unit pairs in the largest component of $\mathbf{A}$. A component is a connected subnetwork such that all units in the subnetwork have infinite path distance to non-members of the subnetwork. This is identical to the bandwidth selection in leunggnn, which verifies the high-level assumptions required to characterize the asymptotic properties of $\hat{\sigma}_n^{ADT} $ under the bandwidth choice (7).

Simulation study

In this section, we conduct Monte Carlo simulations to evaluate the finite-sample properties of our proposed estimator. Its performance is compared against conventional regression models designed to account for spillover effects. The data generation process (DGP) is specified as follows: First, an undirected network $\mathbf{A}$ is generated using the Erd\H{o}s--R\'enyi model, setting the connection probability between any two distinct nodes to $0.5$. Simulations are run for network sizes $n \in \{1500, 2000, 2500, 3000\}$. For each node $i$, 15 covariates are generated from a multivariate normal distribution, $\mathcal{N}(\mathbf{0},\mathbf{I_{15\times 15}})$, and partitioned into 10-dimensional observed covariates $\mathbf{X_i}$ and 5-dimensional unobserved confounders $\mathbf{U_i}$. Finally, independent error terms $(e_i, \eta_i, \vartheta_i)$ are drawn from $\mathcal{N}(0, 0.1)$. A critical feature of our DGP is the incorporation of network structures into the generation of key variables. The treatment assignment $(D_i)_{i=1}^n$, the outcomes $(Y_{i0}, Y_{i1})_{i=1}^n$, and double negative controls variables $(W_{i}, Z_{i})_{i=1}^n$ are all functions of not only the ego's covariates ($\mathbf{X_i}, \mathbf{U_i}$) but also their neighbors' covariates. Let $\beta_X$ = $(0.1, 0.5, 0.5^2, 0.5^3, 0.5^4, 0, 0, 0, 0, 0)$ and $\beta_U$ = $(0.1, 0.5, 0.5^2, 0, 0)$, and define

and $\Delta Y_i$ = $Y_{i1}$ - $Y_{i0}$. The above models imply that the ADT is $\tau$, which we set to be $0.5$. Next, the basic settings for the GNN estimation are explained. Specifically, an exposure mapping is first intentionally misspecified as $G_i=\mathbf{1}\{\sum_{j=1}^{n-1}A_{ij}D_j>0\}$ (the correct specification being $G_i=\sum_{j=1}^{n-1}A_{ij}D_j/\sum_{j=1}^{n-1}A_{ij}$), and $\tau(1)$ is then estimated. Our GNNs model utilizes the principal neighborhood aggregation (PNA) architecture from corso2020principal. This architecture enhances expressive power by integrating multiple aggregators, making it more robust in networks where node degrees vary widely. In particular, we employ “mean” and “sum” aggregators to combine information from neighboring nodes. To stabilize the aggregation process with respect to node degree, these aggregators are paired with three degree scalers: “identity”, “amplification”, and “attenuation”. The “amplification” and “attenuation” scalers are based on the logarithm of the node degree, a design that helps prevent the exponential amplification of gradients across successive GNNs layers. All GNNs estimators are set with $L = 1$. The update function within each layer consists of a MLP with a hidden layer dimension of 16, using $\mathbf{ReLU}$ as the activation function. For model training, we employ the Adam variant of stochastic gradient descent for optimization. The learning rate is fixed at 0.01, and the number of training epochs for each nuisance function is 500. All neural network implementations are based on the $\mathbf{torch}$ package in $\mathbf{R}$. Subsequently, we compare the proposed estimator with the following model:

where both $\beta$ and $\gamma$ are 10-dimensional column vectors. This model uses the correct exposure mapping $G_i=\sum_{j=1}^{n-1}A_{ij}D_j/\sum_{j=1}^{n-1}A_{ij}$ and accounts for the covariates of neighbors. Table (ref) compares our proposed estimator ($\tau_{ADT}(1)$) and Model (9) for sample sizes $n = 1500$ to $3000$. Despite a deliberately misspecified exposure mapping, our estimator demonstrates consistency: as $n$ increases, both bias and RMSE decrease (RMSE drops from 0.3240 to 0.2921), indicating robustness to unmeasured confounding. In contrast, Model (9), even with the correct exposure mapping, exhibits persistent bias ($\approx 0.37$) and stable RMSE ($\approx 0.373$) across all sample sizes. This lack of convergence confirms its inability to address omitted variable bias from $\mathbf{U_i}$, yielding precise but biased estimates.

In summary, the simulation results in Table (ref) provide evidence for the superiority of our proposed method. Even when subject to a misspecified exposure mapping, our GNN-based doubly robust estimator yields significantly lower bias and RMSE compared to the conventional regression-based approach.

Empirical application

In this section, we apply our methodology to estimate the causal effect of China's 2012 Green Credit Policy on corporate green innovation huang2023green,li2024impact. Crucially, this evaluation necessitates explicitly accounting for unmeasured network confounding, as corporate responses are often driven by latent factors?such as informal political connections and executive networks?that influence both regulatory treatment and innovation. Relying solely on observables would likely yield biased estimates. Moreover, given the complex, non-linear nature of knowledge spillovers, our framework integrating double negative controls with GNNs is uniquely suited to robustly proxy for these latent confounders and flexibly model high-dimensional network interactions. We analyze data from Chinese A-share listed companies in 2011 (pre-policy) and 2013 (post-policy), sourced from the China National Intellectual Property Administration (CNIPA) and the China Stock Market $\&$ Accounting Research (CSMAR) Database. After excluding financial institutions and firms with data irregularities (e.g., ST stocks), the final sample comprises 1,664 observations. The network adjacency matrix is defined as $A_{ij} = 1$ if firms share the same industry or are located in adjacent provinces. The outcome $Y_{it}$ is measured as the natural logarithm of granted green patents plus one. The treatment variable $D_{i}$ is defined as follows:

where the six major heavily polluting industries include thermal power, iron and steel, petrochemicals, cement, non-ferrous metals, and chemicals liu2019green. Based on the treatment status $\mathbf{D}$, let $W$ (the NCO) and $Z$ (the NCE) be the proportion of treated individuals among an individual's first-order neighbors and non-neighbors, respectively. Notably, this choice of a double-negative control is consistent with Assumption 2.3. In addition, we include 10 control variables in the model, including the proportion of institutional investor shareholding (InstHold), duality of CEO and chairman (Duality), proportion of independent directors (IndDir), ratio of owner's equity to market value (EqToMV), asset-liability ratio (Lev), ratio of cash to total assets (CashTA), ratio of total fixed assets to total assets (FixTA), return on total assets (ROA), number of employees (EmpNum), and Tobin's Q (TobinQ). Table (ref) reports the descriptive statistics of the main variables used in the study.

Next, we select different exposure mappings to explore the impact of GCP on enterprises' green innovation in Table (ref). The first mapping $G_{1i}=\mathbf{1}\{\sum_{j=1}^{n-1}A_{ij}D_j>0\} = 1$ indicates that there exists a treated group around the enterprise. The second mapping $G_{2i}=\mathbf{1}\{\sum_{j=1}^{n-1}A_{ij}D_j>1\} = 1$ indicates that there are at least two treated groups around the enterprise. And the third mapping $G_{3i} = \mathbf{1} \left( \sum_{j=1}^{N} A_{ij} D_j > \frac{1}{5} \left[ \frac{1}{N} \sum_{k=1}^{N} \left( \sum_{j=1}^{N} A_{kj} D_j \right) \right] \right)=1$ indicates that the number of treated neighbors of an enterprise is greater than one-fifth of the average number of treated neighbors across all enterprises. The first column, $ADT(G_{1i}=1)$, estimates the direct policy effect for the subset of treated firms that have at least one treated neighbor. The resulting ADT is 0.2662 and is highly statistically significant (P-value = 0.0014). This indicates that the GCP had a significant direct impact on heavily polluting firms, strongly incentivizing or compelling them to pursue green innovation, particularly when they were networked with at least one other polluting firm. The second column, $ADT(G_{2i}=1)$, focuses on treated firms with at least two treated neighbors. The estimated ADT increases in magnitude to 0.3074. Numerically, this suggests that the direct effect of the policy is even stronger for treated firms embedded in denser clusters of their peers. Finally, the third column, $ADT(G_{3i}=1)$, uses a relative exposure measure. It shows the ADT for treated firms whose number of treated neighbors is greater than one-fifth of the network average. The estimate is 0.1405 and remains statistically significant at the 5$\%$ level (P-value = 0.0251). This confirms that the significant direct effect of the GCP on polluting firms is robust, particularly for those located in relatively “dense” clusters of other treated enterprises.

In summary, the results in Table (ref) reveal significant heterogeneity in the direct effect of the GCP. The policy's direct impact on promoting green innovation among polluting firms is robust and significant, and this direct effect appears to be (at least in magnitude) stronger for treated firms that are more deeply exposed to other treated firms within their network.

Conclusion

Estimating causal effects from observational network data is challenging due to unmeasured network confounding. Traditional methods, including standard DID, often fail in these settings. Therefore, we propose a general DID framework integrating DNC and GNNs. We establish an identification strategy founded upon the latent parallel trends assumption, leveraging DNC via confounding bridge functions. This approach yields doubly robust, non-parametric identification of the treatment effects, requiring only one of the two bridge functions to be correctly specified for identification to hold. Subsequently, we develop a doubly robust estimator. This estimator leverages GNNs combined with the GMM to flexibly estimate the high-dimensional bridge functions inherent in the DNC approach. We further establish the estimator's asymptotic normality under conditions of $\psi$-network dependence and ANI, and provide a network HAC variance estimator. Simulation studies confirm the favorable finite-sample performance of our estimator relative to conventional models, even under misspecification of the exposure mapping. Finally, an empirical application evaluating the impact of China's green credit policy demonstrate the method's practical utility.

Conflict of interest

On behalf of all authors, the corresponding author states that there is no conflict of interest.

Funding

This work was supported by the [National Social Science Fund of China], [24BTJ064].

A Proofs

A.1 Proof of Theorem 1

This proof adopts the methodology of miao2018identifying and adapts it to our DID framework under network interference. First, we will prove

Under Assumption 2.1 and 2.2,

where the first equality follows from Assumption 2.1 and the second equality follows from Assumption 2.2. Under Assumption 2.3,

where the equality follows from Assumption 2.3(a). Under Assumption 2.4, $E[\Delta Y_{i}|G_{i}=g, D_{i}=d, \mathbf{U}, \mathbf{X}, \mathbf{A}] = E[h_1( W_{i}, \mathbf{X}, \mathbf{A})|G_{i}=g,D_{i}=0, \mathbf{U}, \mathbf{X},\mathbf{A}]$, we get

Next, we prove that the confounding bridge function is identified by

using Assumption 2.3, we have

Similarly, we have

Then, under Assumption 2.4, we proof

Finally, under Assumption 2.5, we show that the unique solution to equation (A.1) identifies the outcome confounding bridge function $h$ . Suppose there are two functions $h_1(W_{i}, \mathbf{X}, \mathbf{A})$ and $h'_1(W_{i}, \mathbf{X}, \mathbf{A})$ that satisfy equation (A.1). Then,

for all g, $\mathbf{x}$, $\mathbf{A}$, and almost all $z$. Then, under Assumption 2.5, $h_1(W_{i}, g, \mathbf{X}, \mathbf{A})$ = $h'_1(W_{i}, g, \mathbf{X}, \mathbf{A})$ almost surely. Thus, the solution to equation (A.1) identifies the outcome confounding bridge function.

A.2 Proof of Theorem 2

Similar with Theorem 1, we first prove

Under Assumption 2.1 and 2.2,

where the first equality follows from Assumption 2.1, the second equality follows from Assumption 2.2, and the last follows from the law of iterated expectations. Under Assumption 2.3,

where the third equality follows from Assumption 2.3(b), $Z_{i}\perp \Delta Y_{i} \mid D_{i},G_{i},\mathbf{U}, \mathbf{X}, \mathbf{A}$. Under Assumption 2.6, $E[\mathbf{1}_i(g,1)|\mathbf{U}, \mathbf{X}, \mathbf{A}]/ E[\mathbf{1}_i(g,0)| \mathbf{U}, \mathbf{X}, \mathbf{A}] = E[q_1( Z_{i}, \mathbf{X}, \mathbf{A})|G_{i}=g,D_{i}=0, \mathbf{U}, \mathbf{X},\mathbf{A}]$, we get

Next, we prove that the confounding bridge function is identified by

using Assumption 2.3, we have

Similarly, we have

Then, under Assumption 2.6, we proof

Finally, under Assumption 2.7, we show that the unique solution to equation (A.2) identifies the outcome confounding bridge function $h$ . Suppose there are two functions $q_1(Z_{i}, g, \mathbf{X}, \mathbf{A})$ and $q'_1(Z_{i}, g, \mathbf{X}, \mathbf{A})$ that satisfy equation (A.2). Then,

for all g, $\mathbf{x}$, $\mathbf{A}$, and almost all $w$. Then, under Assumption 2.7, $q_1(Z_{i}, \mathbf{X}, \mathbf{A})$ = $q'_1(Z_{i}, \mathbf{X}, \mathbf{A})$ almost surely. Thus, the solution to equation (A.2) identifies the outcome confounding bridge function.

A.3 Proof of Proposition 1

Under $\mathcal{M}_1$ where working model $h^*_1(W_{i}, \mathbf{X}, \mathbf{A})$ is correctly specified, we have $h^*_1(W_{i}, \mathbf{X}, \mathbf{A})$ = $h_1(W_{i}, \mathbf{X}, \mathbf{A})$ for $g$ $\in$ $\mathcal{G}$. We have

where the third equation is obtained by the law of iterated expectations, the fourth equation follows from

and the last equation follows from the equation (5) in Theorem 1. Based on the identification result of Theorem 1, we have $\tau^{dr}_{ADT}(g)$ = $\tau_{ADT}(g)$. Under $\mathcal{M}_2$ where working model $q^*_1(Z_{i}, \mathbf{X}, \mathbf{A})$ is correctly specified, we have $q^*_1(Z_{i}, \mathbf{X}, \mathbf{A})$ = $q_1(W_{i}, \mathbf{X}, \mathbf{A})$ for $g$ $\in$ $\mathcal{G}$. We have

where the third equation is obtained by the law of iterated expectations, the fourth equation follows from

and the fifth equation follows from the equation (6) in Theorem 2. Based on the identification result of Theorem 2, we have $\tau^{dr}_{ADT}(g)$ = $\tau_{ADT}(g)$. Together, we complete the proof.

A.4 Proof of Theorem 3

This proof adopts the methodology of leung2022 and adapts it to our framework. Let $\mathcal{F}_n$ be the $\sigma$-algebra generated by $(\mathbf{Z}, \mathbf{X}, \mathbf{A})$, either $h, h' \in \mathbb{N}$, $f \in \mathcal{L}_h$, and $f' \in \mathcal{L}_{h'}$, or $h = h' = 1$ and $f = f' = i^*$, $s > 0$ and $(H, H') \in P(h, h'; s)$. Define $\xi = f(\mathbf{C_H})$ and $\zeta = f'(\mathbf{C_{H'}})$. Let $\mathbf{D'}, \mathbf{D''}$ each be independent copies of $\mathbf{D}$. Define $\mathbf{D_i^{(s,\xi)}}$ = $(\mathbf{D}_{N_{\mathbf{A}}(i,s)}, \mathbf{D}'_{N^c_{\mathbf{A}}(i,s)})$, $\mathbf{D_i^{(s,\zeta)}}$ = $(\mathbf{D}_{N_{\mathbf{A}}(i,s)}, \mathbf{D}''_{N^c_{\mathbf{A}}(i,s)})$, $\mathbf{D_{-i}^{(s,\xi)}}$ = $(\mathbf{D}_{N_{\mathbf{A}}(i,s)}, \mathbf{D}'_{N^c_{\mathbf{A}}(i,s)})\texttt{\textbackslash}D_i$, $\mathbf{D_{-i}^{(s,\zeta)}}$ = $(\mathbf{D}_{N_{\mathbf{A}}(i,s)}, \mathbf{D}''_{N^c_{\mathbf{A}}(i,s)})\texttt{\textbackslash}D_i$, and

where $\mathbf{1}^{(s,\xi)}_i(g,d)$ = $\mathbf{1}\{G(i, \mathbf{D_{-i}^{(s,\xi)}}, \mathbf{A})=g, D_i=d\}$ for $d$ = 0, 1. The definition of $C_i^{(s,\zeta)}$ is similar to that of $C_i^{(s,\xi)}$. Finally let $\xi^{(s)}$ = $f((C_i^{(s,\xi)} : i \in H))$ and $\zeta^{(s)}$ = $f((C_i^{(s,\zeta)} : i \in H'))$. Based on Assumptions 3.2 and 3.3, $C_i$ is uniformly bounded, so $|Cov(\xi,\zeta|\mathbf{Z}, \mathbf{X}, \mathbf{A})| \le 2||f||_{\infty}||f'||_{\infty}$, so for $s$ $\le$ 2max$\{K,1\}$, we have $|Cov(\xi,\zeta|\mathbf{Z}, \mathbf{X}, \mathbf{A})| \le \Psi_{h,h'} (f ,f')$. Now consider $s$ > 2max$\{K,1\}$, so that $\ell_{\mathbf{A}}(H,H')$ > 2max$\{K,1\}$. By Assumption 3.4, $(C_i^{(\lfloor s/2 \rfloor,\xi)} : i \in H)$ $\perp$ $(C_i^{(\lfloor s/2 \rfloor,\zeta)} : i \in H') | \mathbf{Z}, \mathbf{X}, \mathbf{A}$. Then

where the last line uses the fact that, by Assumption 3.4 (b), $\mathbf{1}^{(s,\xi)}_i(g,d)$ = $\mathbf{1}_i(g,d)$ for $d$ = 0, 1, and by Assumption 3.5, $\max_{i\in N_n}E[|\Delta Y_{i}(\mathbf{D})-\Delta Y_{i}(\mathbf{D}_i^{(\lfloor s/2 \rfloor,\xi)})|\,\big|\, \mathbf{Z}, \mathbf{X}, \mathbf{A}]\le \theta^{ADT}_{n,\lfloor s/2 \rfloor}$.

A.5 Proof of Theorem 4

Decompose

where

For $C_{i1}$, there exist universal constants C > 0,

where the first inequality uses Assumption 3.3, the second inequality uses Lemma B.1, and the last inequality uses Cauchy-Schwarz inequality. The last line is $o_P$(1) by Assumption 3.1(a) and 3.6 (b). For $C_{i2}$, we directly have that $C_{i2}$ = $o_P(1)$ by Assumption 3.1 (c). For $C_{i3}$,

by H$\ddot{\text{o}}$lder's inequality, Assumption 3.1 (b) and Assumption 3.3. The proof of $C_{i4}$ and $C_{i5}$ are similar with $C_{i1}$. For $C_{i4}$, we have

by Assumption 3.3, Assumption 3.6 (b), Lemma B.1 and Lemma B.2. For $C_{i5}$,

by Assumption 3.1 (a), Assumption 3.3, Assumption 3.6 (b), Lemma B.1 and Lemma B.2. The proof of $C_{i6}$ and $C_{i7}$ are similar with $C_{i3}$. For $C_{i6}$,

by H$\ddot{\text{o}}$lder's inequality, Assumption 3.1 (a), Assumption 3.3 and Lemma B.2. For $C_{i7}$,

by H$\ddot{\text{o}}$lder's inequality, Assumption 3.1 (b), Assumption 3.3 and Lemma B.2. Totally, we have

By Theorem 2, $\{C_i\}_{i=1}^n$ is conditionally $\Psi$-dependence given $(\mathbf{Z}, \mathbf{X}, \mathbf{A})$ with the dependence coefficients $\{\widetilde{\theta}^{ADT}_{n,s}\}_{s \ge 0}$. Then, letting $\tilde{C}^{ADT}_{N_n}$ = $n^{-1/2}\sum_{i\in N_n}\left(C_i-\tau^{mr}_{ADT}(g)\right)/\sigma_n^{ADT}$, the same arguments as in the proofs of Lemmas A.2 and A.3 of Kojevnikov et al. (2021) show that there exists a positive constant C > 0 such that

where $\phi$ denotes the cumulative distribution function of $\mathcal{N}(0,1)$, $m_n$ and $\varepsilon$ are as given in Assumption 3.6. The right-hand side converges to zero by Assumption 3.6, implying that $\tilde{C}^{ADT}_{N_n}$ $\overset{d}{\rightarrow} \mathcal{N}$(0, 1). Thus, we have

A.6 Proof of Theorem 5

This proof is similar with leunggnn. Define

where $\tau_{ADT,i}(g)$ = $C_i-\tau^{dr}_{ADT}(g)$. We first show that |$\hat{\sigma}_n^{ADT}$ - $\tilde{\sigma}_n^{ADT}$| $\overset{p}{\rightarrow} $ 0. We have

by H$\ddot{\text{o}}$lder's inequality. Next, for some constant C > 0,

by Assumptions 3.2, 3.3, 3.7(a) and (d). Then, by Assumptions 3.2, 3.3, 3.7(a) and (b), we have

Next, the proof of Theorem 4 of leunggnn can be applied to show that

The argument follows from substituting $\tilde{\tau}_{ADT,i}(g)$ for $Z_i$ - $\tau_i(t,t')$ in Theorem 4 of leunggnn and our Assumptions 3.7(c)-(e) for Assumptions 7(b)-(d). Finally, we apply Proposition 4.1 of kojevnikov2021limit to show $|\hat{\sigma}_n^{ADT*}-\sigma_n^{ADT}| \overset{p}{\rightarrow}$ 0.

B Lemmas

$\mathbf{Lemma\ B.1.}$ Under Assumptions 3.2-3.5,

for $d$ = 0, 1, a constant 0 < $\varepsilon$ < 1, holds with $\widetilde{\theta}^{ADT}_{n,s}$ = $\theta^{ADT}_{n,\lfloor s/2 \rfloor}\mathbf{1}\{s > 2max\{K,1\}\}+\mathbf{1}\{s \le 2max\{K,1\}\}$ for all $n$ $\in$ $\mathbb{N}$ and $s$ > 0 and

for either $h$, $h'$ $\in$ $\mathbb{N}$, $f$ $\in$ $\mathcal{L}_h$, $f'$ $\in$ $\mathcal{L}_{h'}$, or $h$ = $h'$ =1 and $f$ = $f'$ = $i^*$. $\mathbf{Proof:}$ Similar with the proof of Theorem 1, we have

then, the proof was completed using Corollary A.2 of kojevnikov2021limit. $\mathbf{Lemma\ B.2.}$ Suppose that Assumptions 3.2, 3.4 and 3.6 hold. Then, we have

for all $g$ $\in$ $\mathcal{G}$ and $d$ = 0, 1. $\mathbf{Proof:}$ By Assumptions 3.4 (c), we have $\mathrm{E}[\hat{\mathrm{E}}[\mathbf{1}_i(g,d)]]$ = $\mathrm{E}[\mathbf{1}_i(g,d)]$, and thus it suffices to show that $Var(\hat{\mathrm{E}}[\mathbf{1}_i(g,d)])$ = $O(\frac{1}{\sqrt{n}})$. Observe that

where the last equality follows from Assumption 3.4. By the Cauchy-Schwarz inequality, the second term of the last line is bounded above by $n^{-1}\sum_{s = 1}^{2K}M^{\partial}_{N_n}(s)$ which is $O(n^{-1})$ by Assumption 3.6 (a).

C Average indirect effect

The proof of AIT follows a similar approach to that of ADT. Therefore, in the following, we present the propositions and theorems related to AIT and omit the proof. First, we propose the following identifying assumption for AIT: $\mathbf{Assumption\ C.1.}$ \\ (a) (Latent parallel trends for AIT)

(b) (Negative controls) \\ Negative control outcome (NCO): For all $i$ $\in$ $N_{n}$ and all $j$ $\in$ $\mathcal{E} _{i}$, $W_{j}$ satisfy

Negative control exposure (NCE): For all $i$ $\in$ $N_{n}$ and all $j$ $\in$ $\mathcal{E} _{i}$, $Z_{i}$ satisfy

$\mathbf{Assumption\ C.2.}$ \\ (a) (Outcome confounding bridge function) There exists a function $h_2$($W_{j}$, $D_{i}$, $\mathbf{X}$, $\mathbf{A}$), such that for all $i$ $\in$ $N_{n}$ and all $j$ $\in$ $\mathcal{E} _{i}$,

(b) (Negative control relevance) For any square integrable function $f$ and any $d$, $\mathbf{x}$ and $\mathbf{A}$, if E($f(W_{j}$) $\mid$ $Z_{i}$ =$z$, $D_{i}$ = $0$, $\mathbf{X}$ = $\mathbf{x}$, $\mathbf{A}$) = 0 for almost all $z$, then $f(W_{j}$) = 0 almost surely. Then, we establish non-parametric identification of the AIT under Assumption 2.2, C.1 and C.2. $\mathbf{Theorem\ C.1}$ Under Assumption 2.2, C.1 and C.2, the confounding bridge function is identified as the unique solution to the following equation:

and the AIT is identified by

In a DAG framework similar to Figure 2, given a sufficiently large sample size, $D_{j}$ satisfies the NCO conditions, while the treatments of units located at a distance of at least 2 from individual $j$ $\in \mathcal{E} _{i}$ fulfill the NCE requirements. Second, we givean alternative identification approach and subsequently derive the doubly robust DID estimands for AIT.

$\mathbf{Assumption\ C.3.}$ \\ (a) (Treatment confounding bridge function) There exists a function $q_2$($W_{j}$, $\mathbf{X}$, $\mathbf{A}$), such that for all $i$ $\in$ $N_{n}$ and all $j$ $\in$ $\mathcal{E} _{i}$,

(b) (Negative control relevance) For any square integrable function $f$, $\mathbf{x}$ and $\mathbf{A}$, if E(f($Z_{j}$) $\mid$ $D_{i}$ = $0$, $W_{i}$ =$w$, $\mathbf{X}$ = $\mathbf{x}$, $\mathbf{A}$) = 0 for almost all $w$, then f($Z_{i}$) = 0 almost surely. we establish another identification of the AIT under Assumption 2.2, C.1 and C.3. $\mathbf{Theorem\ C.2}$ Under Assumption 2.2, C.1 and C.3, the confounding bridge function is identified as the unique solution to the following equation:

and the AIT is identified by

Let $h^*_2(W_{j}, \mathbf{X}, \mathbf{A})$ and $q^*_2(Z_{j}, \mathbf{X}, \mathbf{A})$ be arbitrary models for the true, unknown function $h_2(W_{j}, \mathbf{X}, \mathbf{A})$ and $q_2(Z_{j}, \mathbf{X}, \mathbf{A})$ for $i$ $\in$ $N_n$ and $j$ $\in$ $\mathcal{E} _{i}$. In order to describe our proposed doubly robust approach, consider the following two models:

The doubly robust DID estimand can be expressed as:

The following Proposition C.1 shows that our proposed doubly robust estimand recovers the ADT provided that at least one of models $\mathcal{M}_1$ and $\mathcal{M}_2$ is correctly specified. $\mathbf{Proposition\ C.1}$ \\ If at least one of the models $\mathcal{M}_3$ and $\mathcal{M}_4$ is correctly specified, then $\tau^{dr}_{AIT}$ = $\tau_{AIT}$.