Neighborhood Adaptive Estimators for Causal Inference under Network Interference

Neighborhood Adaptive Estimators for Causal Inference under Network Interference

Introduction

Estimation of causal effects has become a staple, if not a requirement, in many applied fields. This breadth comes due to the realization that scientific and societal knowledge requires rigor beyond predictions. This boon in applications has strained existing causal inference tools and led to a mass development of theory and methodology that can address modern causal questions. Much of the early work on causal inference concerned itself with the design and analysis of controlled experiments where a plethora of independence assumptions made them plausible and implementable splawa1990application,fisher1935design. This work was then generalized by developing classes of assumptions and methodologies that allowed practitioners to emulate simple experimental designs using observational data, still under various independence assumptions rubin1978bayesian,manski1990nonparametric. However, many modern causal settings have made these assumptions untenable: individuals, corporations, and countries are increasingly connected, with the deployment of individual treatments potentially affecting the outcomes of large groups of individuals, either directly or indirectly.

For example, an individual's probability of changing housing may increase if (i) he or she receives a direct incentive to move or (ii) if his or her nearby households receive a housing voucher to move sobel2006randomized. Similarly, the probability that an individual will be infected by a virus can decrease by several mechanisms: (i) he or she can receive the vaccine, which provides direct protection against infection, or (ii) his or her friends could be vaccinated, making them less likely to be infected, which in turn reduces the amount of exposure on their friends ross1916application,zivich2021assortativity. In these settings, there is a potential direct treatment effect (from a unit's own treatment status) and a potential indirect treatment effect (due to the treatment status of other units); see cai2015social,hu2021average for other examples. Importantly, the possibility of this indirect treatment effect, dubbed “interference” among units is a violation of the standard Stable Unit Treatment Value Assumption rubin1978bayesian,rubin1980randomization and Individualized Treatment Response manski1990nonparametric assumption. Understanding what estimands are plausible in these settings, what identification assumptions are needed, and what estimation strategies can be employed has been the central focus of recent work. This is particularly critical in the context of observational studies, where confounding can further complicate the identification and estimation processes.

In order to identify and estimate the treatment effects under interference, different types of estimands and causal assumptions have been considered toulis2013estimation,aronow2017estimating,sussman2017elements,jagadeesan2020designs. For example, hudgens2008toward propose a partial interference model in which units are split into groups where interference exists within each group but not outside the group. This assumption offers relatively well-understood dependence structures between one's potential outcomes and the treatments of others tchetgen2012causal. The constant treatment response (CTR) assumption of manski2013identification defines a class of exposure mappings which can be seen as a representation of the whole treatment vector and any other covariates. The exposure mapping is a general framework to capture complex structures that can influence interference. This has recently been specialized for the setting of network interference aronow2017estimating,savje2017average,savje2021causal.

Although it is widely recognized that interference and other “network effects” can depend on various types of network structures cartwright1956structural,awan2020almost, most of the existing literature focuses on neighborhood interference problems karwa2018systematic, sussman2017elements. In these settings, interference is restricted to those who share edges and, importantly, the strength or type of interference does not depend on whether those with direct edges have connections between them. In many applied settings, this interference characterization may not be sufficiently rich to fully capture the variability in potential outcomes as a function of the full treatment vector savje2021causal. We consider a broad generalization of the neighborhood interference approach by studying interference that can occur within a distance $m_i$ of the unit $i$ in a graph; importantly, this interference can depend on the full subgraph of radius $m_i$, and not just on edge counts or density measures. When $m_i=1$ for every unit, this setting recovers the neighborhood interference assumption of many papers in the literature, e.g. sussman2017elements, forastiere2021identification. However, we allow the value of $m_i$ to be unknown, to vary between units, to depend on the treatment assignment configurations, and to possibly be large.

One of the main limitations in estimating causal effects in the presence of interference is the large amount of data that is needed to estimate treatment effects based on flexible exposure mappings. For example, the exposure mapping “(is treated?, has at least one treated friend?)” maps to four possible counterfactuals, while the naive generalization to “(is treated?, has $k$ treated friends?)” maps to $2 d$ possible counterfactuals where $d$ is the degree of an individual. Even this simple generalization can lead to positivity violations if care is not taken in specifying how the potential outcomes interact with exposure mappings. With that in mind, our aim is to develop adaptive estimators for interference that balance the availability of data and the interference patterns themselves. That is, we will consider unit-specific values of $m_i$ that control the complexity of individual potential outcomes.

This approach is in the spirit of the Constant Treatment Response assumption of manski2013identification which suggests that different individuals may have different global treatment patterns that correspond to the same counterfactual values. Thus, our work aims to contribute to the somewhat fundamental question of obtaining an adaptive restriction of the potential outcomes that allows for a more efficient estimation of the estimands of interest. In our case, we search for these restrictions based on available network information, treatment assignment configurations, and the induced bias and variance trade-off. The network provides a natural hierarchy to search over, while the treatment assignments (the postulated source of interference) allow us to match patterns across the network. From an estimation perspective, since we are searching across possible patterns, this approach is connected with model selection and feature engineering, where regularization tools play a key role. As discussed in the following, we construct a form of regularization based on synthetic treatments that replaces the role of sample splitting, which is harder in network settings. Moreover, our approach incorporates possible misspecifications on model selection and controls their impact to obtain adaptive estimators.

Our setting follows the additive effect under the Neighborhood Interference Assumption (ANIA) framework in sussman2017elements. sussman2017elements showed that the potential outcomes satisfy additivity of main effects if and only if the potential outcomes can be parameterized through a partially linear model. Although we have different results for the estimation of the average direct treatment effect on the treated, our main result pertains to the case where there is selection on the treatment assignment but those are independent across units given observables.

Recent development beyond neighborhood interference

Our proposed approach defines an explicit hierarchy of interference patterns in a network, allowing an estimator to trade off bias and variance of the estimation of counterfactual outcomes. We note that the general notion that interference can be hierarchical and reach deep through a social network has previously been mentioned, albeit in passing and with the requirement that interference is the same for each individual within the network puelz2019graph,eckles2016design,jagadeesan2020designs. Although only briefly mentioned, the generalization of these approaches is natural: if friends can lead to interference, then why not friends of friends? However, on closer inspection, it is not immediately clear how this would be implemented. For example, in jagadeesan2020designs this would require specifying the smoothness of the interference function that would presumably depend on the depth of interference. This may or may not lead to a specification with a computationally tractable design. Similar issues with computation are likely when you need to calculate the statistic in tests for deeper than one neighborhood network interference.

Recently, leung2022causal proposed an approximate neighborhood interference (ANI) assumption that states that interference decays when the distance between treatments and the ego increases. When the distance goes to infinity, the interference decays to zero. In this approach exposure mappings are used to define interesting estimands that summarize treatment and spillover effects. An explicit connection between this assumption and a type of generative network dependence called $\psi$-dependence leads to the result that the Horvitz$ $-Thompson estimator for the direct effect is asymptotically normal kojevnikov2021limit. This approximation is very much in the spirit of our work.

Another important development in the literature concerns misspecified exposure mappings. While knowing the correct exposure mapping provides a natural recipe for asymptotically high-quality estimators of many estimands, misspecification of the exposure mapping leads to inappropriate bias in estimation. savje2021causal characterizes this bias and provides conditions under which consistency of estimators can be achieved (essentially a rate of decay of the misspecification of the exposure mapping). Although our work is not guaranteed to be robust to all misspecification, by choosing a rich class of interference patterns, we hope to alleviate many of these concerns.

In the same spirit, yuan2021causal propose a methodology to detect the exposure conditions for a given dataset by precomputing a large class of causal network motifs for each unit (i.e., all possible treated neighborhood graphs for each unit). They then use regression trees to identify important causal network motifs and thus limit the classes of exposures that they need to study. This is similar in spirit to the approach in awan2020almost, which selects relevant causal motifs to adjust for interference in an experiment. Both of these approaches become computationally challenging as the degree of nodes and the depth of interference increase even to modest sizes.

An interesting recent paper leung2022unconfoundedness proposes the use of graphic neural networks (GNN) to learn the relevant representation of interference in a very general setting with possible joint determination of treatment. The paper establishes both positive and negative identification results in general settings. The work also derives asymptotic normality results under appropriate assumptions on the rates of convergence of the GNN predictions. A key related question in that paper is the determination of the relevant depth of interference which directly relates with the number of layers in the GNN. leung2022unconfoundedness highlights the need to impose additional structure to strike some balance between applicability, interpretability, and tractability. Our setting attempts one such balance by imposing assumptions of additivity of the direct treatment effect while allowing nonparametric estimates for the interference function where the relevant depth of interference can be node specific and dependent on the treatment assignments as well.

Our contribution

Our proposed approach differs from the previous literature by recognizing that not all individuals experience interference in identical ways. We codify this by allowing the interference experienced by each individual to depend on a potentially different subgraph of the larger social network. The following examples illustrate the need for different interference maps for different individuals within the same network.

Our broad goal in this work is to study the direct effect of a treatment on individuals while accounting for any potential and variable interference. Under common additivity, positivity and unconfoundedness assumptions and using our new generalization of neighborhood interference (described in detail in Section (ref)), we derive an augmented inverse probability weighting (AIPW) estimator and an outcome regression (OR) estimator for the average direct effect on the treated. We demonstrate that these estimators are asymptotically well behaved when the types of interference are known, and more importantly, even when those need to be adaptively selected within the estimator.

Our analysis of the proposed interference estimator builds upon ideas from the variable length Markov chain literature for discrete alphabet time series; see, e.g., buhlmann1999variable,ferrari2003estimation,garivier2011context,o2012adaptive, belloni2017approximate. In that setting, the order of the Markov process is unknown and can depend on the specific past history. In our setting, for each node, we will learn the number of hops\footnote{Here a $m$-hop neighborhood of node $i$ includes the induced subgraph of all nodes for whom the shortest path to $i$ is less than or equal to $m$ edges.} needed to obtain a good approximation for the interference which depends on the (sub)network and assigned treatments. We will rely on a Lepski-type estimator (e.g. lepskii1992asymptotically,lepski1997optimal,birge2001alternative,spokoiny2016nonparametric) to select the relevant order $k$ for each individual that defines the relevant (sub)network interference patterns. However, because the choice of $k$ is allowed to also depend on the treatment assignments to nearby nodes, the learning of such treatment-graph patterns is a form of feature engineering that creates challenges. Thus, the proposed approach is an adaptive extension of many existing estimators that assumes a known and fixed $k$-hop neighborhood determines the interference that complements the literature in that direction. Our interference estimator builds on belloni2017approximate which estimates a variable length Markov chain for a stationary time series. In our setting, we establish oracle inequalities for the performance of the interference estimator. This result implies that our interference estimator performs similarly to an (oracle) estimator that optimally balances (the unknown) bias and variance within the class. In well-specified (low-degree) cases, we provide specific rates implied by the oracle inequality.

As mentioned earlier, the estimation of the interference function is an intermediary step in the estimation of the average direct treatment effect (ADTE). The analysis of the outcome regression estimator is a direct consequence of the estimation of the interference function, and it inherits the rates of convergence. On the other hand, the augmented inverse probability weighted estimator leverages the independence between treatment assignments (while allowing heterogeneous propensity scores) to create an orthogonal moment condition which helps diminish the impact of using estimated nuisance parameters (e.g. interference function estimates). However, our interference estimator not only learns the relevant interference values, but also generates the patterns to be tested. Such feature engineering creates an additional dependence that impacts the rate of convergence of the interference estimator. This dependence is problematic in the analysis of the distribution of the AIPW estimator. To mitigate this dependence, we pursue a way to decouple feature engineering from the rest of the estimation. Sample-splitting has been a standard way to decouple the dependence within nuisance parameters, but, in our setting, the network dependence severely complicates the use of sample-splitting since we are allowing for the dependence radius to be unknown. A major contribution of our work is to propose the so-called synthetic treatments to bypass these difficulties. This leads to a modified interference estimator that breaks the dependency between the feature engineering step and the estimation of the interference function. In our analysis, the generation of synthetic treatments acts similarly to sampling-splitting and cross-fitting (used in cross-sectional data) to reduce the dependence across the nuisance function estimators. Our analysis of the AIPW estimator accounts for the feature engineering step of learning the interference patterns and provides conditions under which the estimator is asymptotically unbiased and asymptotically normal. Finally, we discuss conservative and non-conservative estimates for the variance.

Setup

We consider a finite population model with $n$ units connected through a network $G$ that articulates how interference can propagate between units. We denote by $\widetilde Y_i(z) = \widetilde Y_i(z_i,z_{-i})$ the potential outcome of the $i$th unit when $z\in\{0,1\}^n$ is the assigned treatment to the whole network.\footnote{Under the standard SUTVA, which we do not make, one would have $\widetilde Y_i(z)=\widetilde Y_i(z_i,0_{-i})$ for all $z\in\{0,1\}^{n}$.} We observe the network $G$ (interference graph), and the triple of treatment, covariates and outcomes. $(Z_i, X_i, Y_i)_{i=1}$, for each unit $i \in [n]$. The outcomes are related to the potential outcomes by

We will consider the following assumptions.

Assumption (ref) is standard in the literature on network interference toulis2013estimation,sussman2017elements,karwa2018systematic,awan2020almost,jagadeesan2020designs. Conditional on covariate information, Assumption (ref) implies that treatment assignment is independent of potential outcomes, while Assumption (ref) implies that treatment assignments are independent. Both of these assumptions are trivially satisfied in the experimental setting, where the treatment is assigned uniformly at random to units. Assumption (ref) rules out a direct effect of $X$ on the potential outcomes of the untreated units. When $X_i$ is of low dimension, it is practical to stratify the analysis by the values of $X$. In that case Assumption (ref) can be relaxed and ideas developed here would also extend to that (admittedly with larger sample size requirements). Indeed, because we do not need to identify their impacts separately, we would still be able to recover the estimand of interest. Nonetheless, Assumptions (ref), (ref) and (ref) allow designs with a nonconstant propensity score $ e(X_i) := \mathbb{P}(Z_i=1\mid X_i)$ and heterogeneous treatment effects $\tau_i$ that depend on $X_i$.

Assumption (ref) also allows for heterogeneous direct treatment effect and is made to simplify the exposition and technical results but can be relaxed in different ways. For example, if covariates are discrete it is possible to estimate the interference function for each covariate value and our results would carry over with minor adjustments. In Section (ref) we discuss a modified estimator that allows for different specifications that postulate the interference to be a (smooth) function of the covariance {\it given an interference pattern}. Moreover, if the covariates $X_i$'s are independent of the graph $G$, the interference is likely to be approximately well balanced between the treatment and control groups. However, we are motivated by situations where we would expect some relation between the covariate $X_i$ and the graph $G$ (for example, $X_i$ is the degree of unit $i$).

Next, we impose restrictions on the potential outcomes. In particular, how the interference network $G$ defines the interference patterns in the potential outcomes. We will allow the interference to depend on a $m$-hop neighborhood of the labelled graph\footnote{Formally we define $G^z$ as the network $G$ where the $i$th node is labelled as $Z_i$. We let $G_i^z$ as $G^z$ but the $i$th node is labeled as the “ego" that is the center to grow a $m$-hop subgraph.} $G^z$ but the value of $m$ is unknown and potentially dependent on the local network and treatments. Here, a $m$-hop neighborhood of node $i$ includes the induced subgraph of all nodes for whom the shortest path to $i$ is less than or equal to $m$ edges.

Assumption (ref) is a restriction to the potential outcomes in the spirit of the constant treatment response assumption of manski2013identification. Without loss of generality, it postulates that the labeled subgraph centered on $i$ with radius $m_i$ is sufficient to determine potential outcomes. Throughout the paper, we have $$\gamma(G_i^z) = \gamma_0(m_i\textrm{-hop neighborhood from }i\textrm{ in }G^z)$$ where $\gamma_0$ is a known mapping and $m_i=m(G_i^{z})$ is an unknown function. Importantly, we allow $m_i$ to be equal $n$ for every $i\in[n]$\footnote{In the next section we introduce approximation errors of using a smaller value than $m_i$ to estimate the interference function for a pattern $g$.}. As discussed below, combined with the additivity assumption (ref), the interference will be a function of $\gamma(G_i^Z)$, that is, $f_i:=f(\gamma(G_i^Z))$. Although the function $\gamma_0$ might be the identity, in applications $\gamma_0$ might summarize the corresponding subgraph. We will require that if the $m$-hop neighborhoods of two units $i$ and $j$ match, then they also match the number of treated nodes in each $k$ hop neighborhood for $k\leq m$. That is, letting $G_i^Z(k):=k\mbox{-hop neighborhood from $i$ in $G^Z$}$, we have that $\gamma_0(G_i^Z(m))=\gamma_0(G_j^Z(m))$ implies that $\gamma_0(G_i^Z(k))=\gamma_0(G_j^Z(k))$ for $k\leq m$. In what follows, $T_{i,k}$ denotes the number of treated nodes at distance $k$ from $i$.

From an interpretation and algorithmic point of view, it is convenient to represent this interference function as a tree. Each local pattern $\gamma_0(G_i^Z(m))$ (built from $m$-hops from unit $i$) corresponds to a node of the tree that is at depth $m$ from the root. The root corresponds to zero hops (e.g., no interference) and the leaves correspond to a relevant interference pattern for the function. (We will associate an interference estimate to each node of the tree, see Section (ref)). Figure (ref) considers an example where $\gamma_0$ is the sum of treated nodes in the relevant local subnetwork $G_i^Z(m)$. The left panel illustrates an interference function that depends only on its neighbors. Figure (ref) on the right panel illustrates an interference function that depends on the “neighbors of the neighbors" if no (immediate) neighbors were treated.

Thus, this assumption is related to the Neighborhood Interference assumption in the causal inference literature (see, e.g., Assumption 2 in forastiere2021identification). The function $\gamma_0$ defines an interference pattern that is similar to the “exposure mapping” discussed in aronow2017estimating. As we discuss below, when we estimate $m_i$ we are concerned with approximating the corresponding interference, and we will balance a bias and variance trade off of the interference (rather than aiming for a perfect recovery of $m_i$ that would rely on strong separation assumptions).

The combination of Assumptions (ref)-(ref) implies the the following partially linear form for the potential outcomes:

where $\epsilon_i$ is a zero mean term (heteroskedastic), and $f(\gamma(G_i^z))$ is the interference term that absorbs the (constant) baseline, see, e.g., Proposition 1 in awan2020almost and sussman2017elements for similar development. Moreover the individual specific $\tau_i:=\mathbb{E}\left[ \widetilde Y_i(1,0_{-i}) - \widetilde Y_i(0,0_{-i}) \mid X_i \right]$ can be dependent of $X_i$. There is a complication if the shocks $\epsilon_i$'s are allowed to be dependent as differentiating between $f(\cdot)$ and the error term would not necessarily be possible shalizi2011homophily. We will make the assumption that the errors are uncorrelated. First, it is plausible that there is no such excess homophily in the problem (such as when connections are enforced exogenously), and so the no-correlation assumption is reasonable. Furthermore, recent work by mcfowland2021estimating suggests that the structure on the covariance of $\epsilon_i$ still allows us to identify $f(\cdot)$ and the error separately. While knowing the structure a priori is a slight generalization to the no-correlation setting, it substantially complicates the derivations without contributing to the understanding of the proposed method and so it is left to future work.

The main estimand of interest we consider is the average direct treatment effect on the treated, i.e.,

Related to ((ref)) we can define the average direct effect on the treated (ADET) as the average over the treatment effect ((ref)) over the treatment assignments

where the outermost expectation is over the treatment assignments. We provide results for this case as well.

Estimators

In this section, we propose estimators for the average direct treatment effect on the treated that leverage the decomposition ((ref)). These estimators build upon the following high-level ideas: (i) under the partially linear assumption, we can obtain estimates $\hat f_i$ of the interference $f_i:=f(\gamma(G_i^Z))$ based on the control units; (ii) subtracting the estimates of the interference from the observed outcomes, $Y_i-\hat f_i$, up to estimation errors, we recover a setting with no interference for which estimators are available. However, estimation errors here include sampling errors, regularization and model selection biases, and errors from learning the features themselves.

To estimate the interference function, we propose a data-driven choice of a unit-specific value of $\hat m_i$-hops that approximates the relevant neighborhood for the interference function. In this setting, up to logarithmic terms, we establish that the proposed estimate balances bias and variance. \footnote{It is based on a Lepski-type procedure to adaptively select a node specific value $m_i$-hops within the control group.}

Based on interference estimates, we can construct the triple $(Y_i-\hat f_i,Z_i,X_i)$ which can be used to construct different estimators for the direct effect. Under the assumption that the treatment assignments are independent conditional on the covariates, we can leverage the propensity score to obtain another estimator: Letting $\hat e_i$ denote an estimate for the propensity score, the {\it augmented inverse propensity weighted} (AIPW) estimator in our setting is given by

This AIPW estimator inherits several good properties that are well known when SUTVA holds. Nevertheless, the construction of the interference function plays a central role in the analysis of these estimators (in particular the impact of the feature engineering step).

Although our main focus is on the estimator ((ref)), under the additive assumption, despite interference, we can estimate the ADTT without the assumption of independent treatment assignment assumption using an {\it outcome regression} (OR) estimator for ((ref)), namely

We will provide rates of convergence based on oracle inequalities for the interference estimator. As described in the following, we will rely on different interference estimators in ((ref)) and ((ref)). The latter has a slower than root-$n$ rate of convergence which is a consequence of the interference estimator. For the former, we need to modify the interference estimator with the use of synthetic treatments to cope with the dependence in our estimation to achieve asymptotic normality for ((ref)).

Interference Estimator

A key step in our approach is to estimate the interference function, which also needs to identify the relevant interference patterns. Our proposal will estimate a relevant $m$-hop neighborhood associated with each $G_i^z$ that achieves a near-optimal bias and variance tradeoff. More precisely, it is a Lepski-type procedure that aims to balance (estimated) bias and variance when choosing $m$. The estimation will construct a tree representation for the relevant patterns that are relevant to approximate interference as discussed in Figure (ref). The estimation of such tree can be viewed as selecting the appropriate subtree out of the tree that contains all possible interference patterns.

In what follows, we use the following notation. Let $G_i^Z(m):=$ $m$-hop subgraph from $i$ in $G^Z$. In order to define a proper partition of the nodes, we match along all intermediate hops as well, that is, if $\gamma_0(G_i^Z(m)) = \gamma_0(G_j^Z(m))$, it implies that $\gamma_0(G_i^Z(k)) = \gamma_0(G_j^Z(k))$ for $k=1,\ldots,m$. For $g=\gamma_0(G_i^Z(m))$, we denote $par(g) = \gamma_0(G_i^Z(m-1))$ and say that $g'\succeq g$ if $\exists m' \geq m$ is such that $g'=\gamma_0(G_i^Z(m'))$.

\\ \uline{{\bf Algorithm 1. Estimation of interference function}} \\ {\bf Step 0.} Initialize $\lambda > 1$, confidence level $1-\delta$, variance upper bound $\bar{\sigma}^2$, $$ \widehat E_n:= \{ g : \exists i,m \in[n], Z_i=0, g=\gamma_0(G_i^Z(m)) \} \ \ \mbox{and} \ \ \widehat{\mathcal{K}}:=\emptyset $$ {\bf Step 1.} For $g \in \widehat E_n$, let $V_g := \{ i \in [n] : Z_i=0, \exists m \in [n] \mbox{ s.t.} \ g=\gamma_0( G^Z_i(m) )\} $ and define $$ \mbox{$\hat f(g):= \frac{1}{|V_g|}\sum_{i\in V_g} Y_i$} $$ {\bf Step 2.} For each $g \in \widehat E_n$, set $\widehat{\mathcal{K}} \leftarrow \widehat{\mathcal{K}} \cup \{g\} $ if $ \exists g',g'' \in \widehat E_n : g'\succeq g, g'' \succeq par(g) $ s.t. $$ \mbox{$ |\hat f(g') - \hat f(g'')| > \lambda \bar\sigma\sqrt{2\log(2n^2/\delta)}\left(\frac{1}{\sqrt{|V_{g'}|}} + \frac{1}{\sqrt{|V_{g''}|}}\right)$}$$ {\bf Step 3.} To estimate $f(\gamma(G_i^Z))$, $i\in[n]$, we define $$

$$ \\

\uline{ } \\

Estimation of interference is challenging as we need to estimate the relevant $m$-hop neighborhood that provides a good approximation for the interference function. Algorithm 1 is based on the idea that, in the control group, the interference function is the signal that we are trying to recover. In its Step 0 the set $\widehat E_n$ denotes the set of different configurations of relevant patterns in the given labeled graph. Because we have a nested subset of graphs as we increase $m$ in $G_i^Z(m)$, we note that the set $\widehat E_n$ can be represented as a tree starting at the root with $m=0$. The set $\widehat{\mathcal{K}}$ is a subtree of $\widehat{E}_n$ that will collect the patterns selected to estimate the interference function.

For each pattern $g \in \widehat E_n$, Step 1 defines an estimate of the interference $\hat f(g)$.\footnote{Note that this is generalizable for more general estimands beyond conditional mean and with covariates.} Step 2 of Algorithm 1 balances the estimates of bias and variance. Informally, a pattern $g$ in $\widehat{E}_n$ should be selected if there are two deeper patterns that have sufficiently different interference estimates relative to their (square root) variance. That is, a pattern $g$ will be selected into $\widehat{\mathcal{K}}$ only if it is needed to keep the bias small relative to an estimate of the variance. The parameter $\lambda$ is a regularization parameter in the spirit of shrinkage estimators, e.g., lasso bickel2009simultaneous. Finally, for each $i$, Step 3 defines the relevant value of hops $\hat m_i$ as the largest value of $m$ that is still in the estimated tree $\widehat{K}$. In the next section, we provide theoretical guarantees on the performance of the estimator. We note that unless one is willing to impose strong separation assumptions, it is unlikely that one can estimate $m$ uniformly across different data generating processes. However, the misspecification cannot be much larger due to the selection criterion in Algorithm 1.

There are several variants of Algorithm 1 that relax some of the requirements. For example, it is possible to use a pattern-specific estimate of the variance $\hat\sigma_g^2$ (instead of assuming a known $\bar \sigma^2$) and our results (e.g., Theorem (ref)) would hold with small modifications. Moreover, the proposed Step 2 that balances bias and variance relies on concentration inequalities (i.e. through the factor $\sqrt{2\log(2n^2/\delta)}$) which can be conservative. It is also possible to rely on bootstrap approximations to take into account correlations to reduce the penalty levels and still achieve the same theoretical guarantees, see, e.g., belloni2018high. Finally, in Section (ref), we discuss different ways to incorporate covariates which would allow us to relax Assumption (ref).

Decoupling feature engineering based on Synthetic Treatments

Here, we propose a variant that changes the definition of the initial tree $E_n$ that will be critical to establishing distributional limits for the AIPW estimator. To see why, we note that Algorithm 1 can be viewed through the lens of feature engineering. In fact, it uses the data to define the patterns to construct the relevant feature (that is, the appropriate $k$ hop neighborhood of the labeled graph). Although such dependence does not affect the rate of convergence of the interference estimator, it is critical for the distributional results of the AIPW estimator.

The next algorithm will be useful when we want to decouple feature generation from treatment assignments. We will generate a different initial set of patterns to replace $\widehat E_n$ as follows.

\\ \uline{{\bf Algorithm 2. Estimation of interference function with synthetic treatments}} \\ {\bf Step 1.} Consider the original (unlabeled) interference graph $G$.\\ {\bf Step 2.} Assign (new) labels $\widetilde Z^o \in \{0,1\}^n$ independently from $Y$ and $Z$.\\ {\bf Step 3.} Define $ \displaystyle \widetilde E_n^o := \{ g : \exists i, m \in [n], g=\gamma_0(G_i^{\widetilde{Z}^o}(m))\} $ \\ \uline{ {\bf Step 4.} Run Steps 1, 2 and 3 from Algorithm 1 with $ \widetilde E_n^o$ instead of $\widehat E_n$. }\\ \\ This decoupling step reduces overfitting in a similar spirit to sample splitting. Nevertheless, the use of synthetic treatments avoids the difficulties that sample splitting would create due to graph dependence across units. The proposal exploits that we only need to create (synthetic) treatments to achieve independence between pruning and feature generation. We acknowledge that, in principle, this decoupling might miss some potential feature that appeared in the true data. However, for a pattern to be informative, it needs to be seen across multiple nodes, and those are likely to appear if we could resample the data. In practice, the label $\widetilde{Z}^0_i$ can be generated from independent Bernoulli random variables with parameters $\hat{e}(X_i)$ for $i \in [n]$. It is also possible to use a small subset of nodes to obtain a rough estimate of the propensity score.\footnote{Under Assumption (ref)(ii) our main results would also hold if we use a small number of observations $q$, $q \ll n^{1/2}/\{(1+D_n) \log n\}$, to learn a rough estimate for $e(X)=\mathbb{P}(Z=1\mid X)$.}

Theoretical Results

In this section, we establish theoretical guarantees associated with the two causal estimators ((ref)) and ((ref)). We also establish the properties of the interference estimator which is key for our proposal.

Our setting follows the Additive effect under Neighborhood Interference Assumption (ANIA) framework in sussman2017elements as the potential outcomes satisfy additivity of main effects if and only if the potential outcomes can be parameterized through a partially linear model. The next assumption is a summary of assumptions (ref), (ref) and (ref) (presented to facilitate readability) with additional moment conditions on the error term.

Assumption (ref) imposes (conditional) independence on the error terms and subgaussianity. It is a technical assumption that can be relaxed in various ways, but it would take the focus out of the main issues the work is focusing on (e.g., finite moments and limited dependence across units). Assumption (ref) also imposes that the individual treatment effects and the interference function are uniformly bounded.

Next, we define the maximum approximation error associated with the $m$-hop subgraphs.

The approximation is defined conditional on a specific vector of treatment assignment $z$. This term provides a quantification of the misspecification one obtains by restricting the model at a specific number of hops for a specific pattern. We note that the set $\bar V_g$ includes both controls and treated units while its subset $V_g$ (used in Algorithm 1) includes only control units. The approximation error $r^Z(g)$ will be used to establish an adaptive bound for our interference estimator. A common assumption in the literature is that controlling for the $1$-hop neighborhood suffices for no approximation errors. There are a few exceptions in the literature that allow for misspecification. The closest to our setting is leung2022causal which considers an approximation $m \to \infty$ so that $r(g_m)\to 0$ and derives asymptotic normality for the inverse propensity weighted estimator under appropriate dependence decay that covers several cases of interest. In contrast, we exploit the heterogeneous radius $m_i=m(G_i^z)$ across different units that depends on the network and the treatment assignments to control the approximation error. This is of interest in obtaining more efficient estimates of the interference. However, this approach creates a model selection problem as we try to balance bias and variance.

Oracle Inequality for Interference Function Estimator

In this section, we establish an oracle inequality for the estimator of the interference function defined by Algorithm 1. The interference function is defined on all labeled graphs and our objective is to obtain a uniform guarantee across all values $f_i := f(\gamma(G_i^Z))$, $i\in [n]$. We obtain an oracle inequality that shows the proposed estimator $\hat f_i$ that is associated with a specific $\hat m_i$-hop subgraph, achieves a similar performance to the best (unknown) choice of $m$ that optimally balances the unobserved bias (measured by $r(g_m)$) and the variance of $\hat f_{g_m}$.

Theorem (ref) establishes a performance guarantee for the estimation of interference in all units. The interference estimate in Algorithm 1 is based only on the control units, thus the effective sample size is $|V_{g_k}|$ (instead of $|\bar V_{g_k}|$). Importantly, the bound holds uniformly for all units. Crucially, the analysis accounts for the fact that $G_i^Z$ is random (Z dependent), which in principle could lead to $n^22^n$ different possible $G_i^z(k)$, $i,k\in [n]$. The result in Theorem (ref) is adaptive in the sense that for patterns $g_k$ that are observed more often, faster rates are obtained, adapting the choice of $k$ for each $i\in [n]$. To provide explicit rates of convergence, we discuss two examples based on two-dimensional regular lattices.

We close this section with an analogous result for the interference estimator based on synthetic treatments proposed in Algorithm 2. The following oracle inequality is a direct corollary of Theorem (ref).

Corollary (ref) has the patterns generated with $\widetilde{E}_n^o$ instead of $\widehat{E}_n$ of Algorithm 1. By resampling the treatment assignments, Algorithm 2 decouples the generation of features from the pruning of the tree. The pruning steps (i.e., Steps 1 and 2) are based on the original assignments and outcomes (since we do not have synthetic outcomes). Corollary (ref) does not provide improvements in the oracle inequality of Theorem (ref). The advantage of using the estimator based on Algorithm 2 is to reduce dependence within the error and feature engineering that appear in some estimators of the ADTT.

Distributional results for the AIPW estimator

Next, we proceed to construct guarantees for the AIPW estimator for the ADTT. Recall that we apply Algorithm 2 to construct the initial set of patterns for Algorithm 1 (that is, to replace Step 0). We use the resulting estimator of the interference function in ((ref)). Our analysis will leverage the independence of the assignment of treatment between units. Next, we state our conditions for the estimation of nuisance functions and their estimability.

Assumption (ref)(i) places mild conditions on the propensity score estimator. Because treatments are assumed to be assigned independently conditioned on observables, the estimation of the propensity score is standard. Assumption (ref)(ii) is a condition on the randomization of the assignments relative to the interference patterns. It assumes the existence of (an unknown and possibly random oracle) $\mathcal{K}^*$ within $\widetilde E_n^o$ that produces a good approximation error and the associated estimator $\hat f^*$ has a minimum convergence rate to estimate the interference function (that is, slightly faster than $n^{-1/4}$). The quantity $D_n$ bounds the number of egos whose neighborhoods are affected by a single node (i.e., bounds the number of changes in $\{\hat f_j\}_{j\in [n]}$ if one treatment is changed). Although it is plausible for $D_n$ to be bounded (e.g., when $m_i^*$ is uniformly bounded and nodes have bounded degree), it is allowed to grow with the sample size. The results are with high probability on the realization of the treatment assignments ($Z$ and $\widetilde Z$). Thus our condition (and analysis) allows for $\mathcal{K}^*$ to depend on the treatment assignments.

We show how these conditions can be used with the orthogonal moment condition implied by the definition of ((ref)). The analysis has new features to deal with the dependence due to interference, and the feature engineering of the proposed approach. We note that as in other context for high dimensional causal inference belloni2014inference,chernozhukov2017double, sample-splitting can be quite helpful in reducing overfitting. However, attempting a sample splitting approach is more complicated due to the network dependence which is also being estimated. We bypass this issue by proposing Algorithm 2 that helps achieve independence from the feature engineering process and the interference estimation.

Next we establish the main result of this section for the estimation of the ADTT. The asymptotic normality result is obtained by relating our proposed AIPW estimator $\hat\tau^{AIPW}$ in ((ref)) to the (unfeasible) estimator $\hat\tau^{AIPW*}$ that knows the true interference function and the propensity score function, i.e., defined as in ((ref)) with $\hat f_i=f_i$ and $\hat e_i = e_i$, $i\in [n]$.

Theorem (ref) establishes the first-order equivalence between the proposed estimator $\hat\tau^{AIPW}$ and the infeasible AIPW estimator that knows the true interference function. In contrast to the outcome regression estimator, the AIPW estimator explicitly exploits the independence of the assignment of treatment between units. This allows us to obtain the so-called Neyman orthogonality and reduce the first-order bias. The analysis further takes into account model selection and feature engineering. This relies on the derivation of sharp bounds conditioned on observing the neighborhood patterns. This is needed to enhance the adaptivity of the approach that does not require the convergence to a deterministic set of patterns. The variance $\Sigma_n$ is needed for the construction of confidence regions. Corollary (ref) provides guidance on its estimation.

Without additional assumptions on the shocks and on the individual treatment effects we can still construct conservative confidence regions based on $\widetilde{\Sigma}_n$. Under the condition that the individual treatment effects can be approximated by an estimable function of the covariates, we can obtain consistent estimates for $\Sigma_n$ and obtain asymptotically correct coverage. In practice, we recommend to pursue the estimation of $\tau_i=\tau(X_i)=\mathbb{E}[Y_i-f_i\mid X_i]$ using the “data” $\{(Y_i-\hat f_i, X_i)\}_{\{i \in [n]:Z_i=1\}}$ by a predictive machine learning model, e.g., random forest. Since only a slow rate of convergence is needed, the condition on the estimators $\hat \tau_i$ in Corollary (ref) is plausible despite of the use of (consistent) estimates $\hat f_i$ of the interference function.

Next we consider randomization over the treatment which alters the estimand and the construction of the confidence interval. It follows that the proof of Theorem (ref) can be used to establish an asymptotic normality result considering the treatment assignment as random.

Proposition (ref) establishes the limiting distribution with the new centering and a correction on the standard error due to the randomness of the treatment. The next proposition discusses the estimation of the relevant variance.

Proposition (ref) shows that $\widetilde{\Sigma}_n^C$ can be used as a conservative estimate of the variance even if we consider the randomness of the assignment of treatment. Moreover, under consistent estimates of the individual treatment effects and propensity scores, we can construct a consistent estimate for $\Sigma_n^D$.

Rates for the Outcome Regression Estimator

Next, we consider the outcome regression (OR) estimator (ref) for ((ref)). This estimator can be consistent without the assumption of independent treatment assignment assumption. For the OR estimator, we apply Algorithm 1 and use the resulting estimator of the interference function in ((ref)). The next result leverages Theorem (ref) to provide a guarantee for $\hat \tau^{OR}$.

Theorem (ref) is based on the oracle inequality of Theorem (ref). A consequence of the proof of Theorem (ref) is that (for $\lambda = 2$) with probability $1-3\delta$ $$|\hat{\tau}^{OR}-\tau | \leq C\max_{i\in [n]} \min_{k \in [n], g_k=\gamma_0(G_i^Z(k))} r^Z(g_k)+\bar \sigma \sqrt{\frac{\log(n/\delta)}{|V_{g_k}|}} $$ which has the maximum rate for the interference function (moreover, note that the last term in Theorem (ref) cannot determine the rate). However, Theorem (ref) averages the adaptive rates across units. Thus, it allows some units to have a small misspecification while other units have a higher level of misspecification. The estimator $\hat\tau^{OR}$ is in general non-regular due to the regularization bias introduced by selecting the relevant $\hat m_i$-hop by balancing (an estimate of) bias and variance. The performance guarantees of $\hat \tau^{OR}$ are robust with respect to the joint dependence on the assignment of treatment between units.

Experimental Performance of the AIPW Estimator

We evaluated the performance of the estimator \(\hat{\tau}^{AIPW}\) under varying graph structures and interference scenarios. Each simulation is based on a fixed sample size of \(n = 2000\) and are repeated \(M = 10,000\) times. Two graph models are considered: Erd\H{o}s R\'{e}nyi (ER) model with independent edges, each having probability \(p = \frac{4}{n}\) and Barabási\-Albert (BA) model which follows a preferential attachment mechanism (with parameter \(\alpha = 1/2\)) where new nodes are more likely to connect to high degree existing nodes. The BA model produces sparser graphs with heavier tailed degree distributions than the ER model.

We consider a single confounder \(X_i = s_i(1)\), where \(s_i(1)\) is the size of the immediate neighborhood of individual \(i\). The treatment assignment mechanism is \(P(Z_i = 1 \mid X_i) = \frac{X_i}{C_0}\), where \(C_0 = (3/2) \max_{i \in [n]} s_i(1)\). This setup ensures that treatment probabilities are proportional to the size of an individual’s neighborhood. The direct treatment effect is $\tau(X_i) = X_i/2$ (meaning that there is greater heterogeneity in the direct treatment effect for BA graphs due to their degree distributions) and we consider two interference functions throughout the simulation:

which decays exponentially on the differences between successive neighborhood sizes, and

where interference is larger the closer unit $i$ is to a treated unit. Lastly, the additive noise in the simulations is \(\epsilon \overset{i.i.d.}{\sim} N(0, 1.5^2)\) and we set the maximum depth of $\widehat{E}_n$ as \(m_{max}\).

Eight scenarios are constructed to examine the interplay between graph structures and interference types. The combinations are labeled via (Graph Type, Interference Function, $m_{max}$), so for example (ER, exp$(n)$, $n$) means the underlying graph is an ER graph, interference has the exponential decay form of depth $n$, and $m_{max} = n$. When $m_{max} = n$ or is greater than the depth of the exponential decay, we apply the pruning rule from in Algorithm $1$. We let $\lambda = 1.01$ and use the synthetic treatments of Algorithm 2.

We report bias, as well as standard errors and coverage for both $\tau$ and $\tau_{ADET}$ in Table (ref). Standard errors for \(\hat{\tau}^{AIPW}(\cdot)\) are constructed using three formulations, reflecting varying degrees of randomness and levels of conservativeness: (1) \(\hat{se} = ( \widehat{\Sigma}_n/\sum_{i=1}^n Z_i)^{1/2}\), where \(\widehat{\Sigma}_n\) is a plug-in estimator of \(\Sigma_n\) from Theorem (ref), as detailed in Corollary (ref); (2) \(\hat{se}_D := (\hat{\Sigma}^D_n/n)^{1/2}\), accounts for the randomness arising from \(\epsilon\) and \(Z\) by employing the plug-in estimator \(\hat{\Sigma}^D_n\), as defined in Proposition (ref), to estimate \({\Sigma}^D_n\) from Proposition (ref); (3) \(\hat{se}_C := (\widetilde{\Sigma}_n/\sum_{i=1}^n Z_i)^{1/2}\), incorporates \(\widetilde{\Sigma}_n\) from Corollary (ref) and Proposition (ref) to construct a conservative estimate of the uncertainty for both \(\tau\) and \(\tau_{ADET}\). Throughout the simulation, \(\hat{se}\) and \(\hat{se}_D\) rely on regression estimates of the individual treatment effects.

The simulations empirically demonstrate the relationship between the standard errors, with \(\hat{se} < \hat{se}_C\) and \(\hat{se}_D < \hat{se}_C\). Throughout the simulations, the greater heterogeneity in the degrees of the BA graphs leads to greater heterogeneity in the individual treatment effects, which in turn makes $\hat{se}_C$ be a substantially more conservative estimator (this is seen through the substantial overcoverage for those graphs). As established in the theory, $\hat{se}$ is the correct standard error for $\tau$, as demonstrated in proper nominal coverage across simulation settings, and similarly $\hat{se}_D$ that incorporates the randomness in $Z$ is the correct standard error for $\tau_{ADET}$ producing proper nominal coverage across those simulations.

Real data application

We reanalize the results of the field experiment from aronow2017estimating and paluck2016changing, which investigated the effects of anti-conflict interventions on school conflict climates. The experiment involved 28 schools and in each school, 15% of students (capped at 64) were identified as eligible for a bi-weekly anti-conflict training intervention based on gender and grade. Half of eligible students were randomly assigned to this intervention. At the end of the school year, various anti-conflict norms were measured, including the self-reported wearing of orange wristbands, which were rewards for anti-conflict behaviors. Previous work identified this measure as particularly salient for assessing the intervention's efficacy. Social networks were constructed based on students’ nominations of up to 10 peers they spent time with forming an undirected graph where an edge exists if either student nominated the other.

We utilize several estimators, each accompanied by valid variance estimators. For \(\hat{\tau}^{AIPW}\), propensity scores are estimated using the covariates “grade" and “gender," consistent with the experimental design in paluck2016changing. The patterns we considered are the number of extra treated units equal to $0$, $1$ or larger than $1$ by increasing the radius of unit's neighbors. For example, $(0,1, \geq 2)$ refers to the set of units where the number of treated units at $0$-th distance neighbor is $0$, the number of treated units at $1$-st distance neighbor is equal to $1$, the number of treated units at $2$-nd distance neighbor is larger or equal to $2$. We estimate $\hat{\tau}^{AIPW}$ under penalty parameter \(\lambda = 1.01\) respectively.\footnote{Since \(\bar{\sigma}\) in Algorithm 1 is unknown, we estimate the variance of each pattern \(g\) by retaining nodes in \(\widehat{E}_n\) with at least five observations. The upper bound in Step 2 of Algorithm 1 is modified to incorporate these variance as: \( |\hat{f}(g^\prime) - \hat{f}(g^{\prime \prime})| \geq \lambda \sqrt{2\log(2n^2/\delta)} \left({\sigma(g^\prime)}/{\sqrt{|V_{g^\prime}|}} + {\sigma(g^{\prime \prime})}/{\sqrt{|V_{g^{\prime \prime}}|}}\right), \) where \(\sigma(g)\) is estimated by sample variance for units with pattern \(g\), defined as \(\sigma_g =({|V_g| - 1})^{-\frac{1}{2}} \cdot \left( \sum_{i \in V_g} (Y_i - \bar{Y}_g)^2\right)^{\frac{1}{2}} \).} The maximum tree depth before pruning is set to \(m_{max} = 4\).

We also consider the H\'{a}jek estimators for the treated population which are compatible with the exposure mapping introduced in aronow2017estimating: $D_i = d_{ab}$ where $a$ is 0 if unit $i$ is not treated and 1 otherwise and $b$ is 0 if unit $i$ does not have any treated neighbors, and 1 otherwise. We let $\hat{\tau}^{HJ}(d_{1j}, d_{0j})$ be the H\'{a}jek estimator contrasting exposure $d_{1j}$ with $d_{0j}$ for $j\in\{0,1\}$. For comparison, we also include \(\hat{\tau}^{NV}\), which does not adjust for confounding or interference. We derive standard error estimators under the two scenarios that the estimand is $\tau$ (written as \(\hat{se}(rs)\)) and $\tau_{ADET}$ (written as \(\hat{se}(rs)\)). When estimating $\tau_{ADET}$ the standard error estimators are constructed by modifying the HAC variance estimators utilized in wang2024designbasedinferencespatialexperiments and applying the law of total variance to account for the joint distribution of treatments and shocks. For \(\hat{\tau}^{AIPW}\), these are given by $\hat{se}(rs) = (\widehat{\Sigma}_n/\sum_{i=1}^n Z_i)^{1/2}$ and $\hat{se}(rst) = ( \widehat{\Sigma}^D_n / n)^{1/2}$.

Table (ref) reports the estimated direct treatment effects using the original exposure mappings and our adaptive approach. We see that failing to adapt to general confounding and interference patterns leads to smalled direct effect estimates. To further explore the underlying structures for the heterogeneity of interference, we investigate \(\widehat{\mathcal{K}}\) under the less conservative penalty parameters \(\lambda = 0.4, 0.7,\) (presented in Appendix (ref)) and \(1.01\). If the correct exposure mapping is defined by \(d_{11}, d_{10}, d_{01},\) and \(d_{00}\), it would be expected that \(\widehat{\mathcal{K}}\) retains only the first layer. However, as shown in Figure (ref), and Figures (ref) and (ref) in Appendix (ref), nodes in deeper layers are consistently retained. This observation suggests the presence of richer structures governing the heterogeneity of interference in this data application.

Additional Discussion and Extensions

There are many possible extensions of the current baseline model that are of interest. In this section, we point out some extensions that are relatively direct to accommodate (possibly increasing the data requirements), while we mention others that would require fundamentally new results.

Throughout the paper we focus on the estimation of the average direct treatment effect on the treated (ADTT) under a simple specification that allows us to clearly communicate the main insights and key contributions (i.e., adaptivity, oracle inequalities, handle feature engineering, etc). Nonetheless the ideas here can be considered in other settings.

For example, the estimation of the average direct treatment effect (ADTE) is a key estimand in causal inference. The ideas proposed in this work extends naturally to this case. In order to pursue this estimand we consider the potential outcomes model described in Assumption (ref) with the additional condition that $\tau_i=\tau(X_i)$. The additional assumption allows us to exploit the idea that $\tau(X_i)=\mathbb{E}[Y_i-f_i\mid X_i]$ when the $i$th node is treated.

Based on Algorithm 1, we obtain estimates of the interference, $\hat f_i, i\in[n]$. Therefore we can estimate the function $\tau(\cdot)$ by estimating a conditional mean model for the data $(Y_i-\hat f_i,X_i)_{\{i\in[n]:Z_i=1\}}$ pertaining to the treated nodes. Letting $\hat\tau_i:=\hat\tau(X_i)$, we have different estimators for ADTE, e.g., $\widehat{\tau}_{ADTE}^{AIPW} = \frac{1}{n} \sum_{i=1}^n \hat\tau(X_i) + Z_i\frac{(Y_i- \hat\tau(X_i) - \hat f_i)}{\hat e_i} - (1-Z_i)\frac{(Y_i - \hat f_i)}{1-\hat e_i}$.

The proposed estimator for $\tau(\cdot)$ relies on an estimated target due to the use of the interference estimator. Despite of the bias that this introduces, the rates of convergence for the interference can be used to establish rates of convergence for the proposed estimator.

The numerical experiments in Section (ref) illustrates a natural extension of our Algorithm 1 to the setting of discrete covariates (i.e. finite support). Importantly, the theoretical development in the paper can be immediately adapted to this setting. The approximation error in Definition (ref) is now defined over patterns of the form $(x,g)$, requiring us to track $V_{x,g}$ throughout. This necessarily has an effect on the effective sample size for estimating $\hat f(x,g)$ in Algorithm 1 but the results of Theorem (ref) (based on Proposition (ref) in the Appendix) follow, mutatis mutandis. The discrete covariate setting is a special case of the following formulation that can handle continuous and other covariates. We note that stratifying on a discrete covariate and then growing the trees in Algorithm 1 has a direct parallel to first growing a single tree based on graph patterns only, and then enforcing a split on the discrete covariate in the leafs. To see the extension to continuous covariates we define $\hat f(g,\cdot) = \arg\min_{h\in\mathcal{H}}\sum_{i\in V_g} (Y_i-h(X_i))^2$ where $\mathcal{H}$ is a class of function that describes the influence of covariates on the interference patterns. If we grow the tree based on $\hat f(g,\cdot)$ that corresponds to splitting on the $X$'s in the leafs. Unlike the discrete case, however, this more general formulation has a more profound effect on the rates of convergence which impacts the pruning of the tree (i.e. the thresholds that are estimated).

We close by noting that we view the current work as an initial contribution to blend machine learning tools within causal inference under network interference for observational studies. Although outside the scope of the current paper, there are several different variations and extensions. Among the extensions that would require fundamentally new results, we would include to handle unobserved graphs, relax the partial linear setting, and to provide a data-driven estimation for the function of patterns $\gamma$. We view the development of tools to allow proper inference under these alternative models as interesting directions for future research.