Causal Inference Under Approximate Neighborhood Interference

Causal Inference Under Approximate Neighborhood Interference

Introduction

\onehalfspacing

Randomized experiments in settings with network interference have seen increasing use in economics and the social sciences.\footnote{E.g.\ bandiera2009social, bond201261, bursztyn_understanding_2014, miguel2004worms, paluck2016changing.} This paper develops methods for causal inference in experiments under relatively weak restrictions on interference. We consider a finite population of $n$ units connected through a network $\bm{A}$. Let $Y_i(\bm{d})$ denote the potential outcome of unit $i$ under the counterfactual that the network is assigned treatment vector $\bm{d} = (d_i)_{i=1}^n \in \{0,1\}^n$. The dependence of $Y_i(\bm{d})$ on the entire vector of assignments allows for interference or “spillovers,” in contrast to the standard potential outcomes model. Interest centers on “exposure effects,” defined below, which summarize how outcomes change in response to manipulations of $\bm{d}$.

The main inferential challenge is that, with a single network, the econometrician observes only one realization of the treatment assignment vector $\bm{D} = (D_i)_{i=1}^n$, where $D_i \in \{0,1\}$ denotes unit $i$'s realized assignment. Identification of exposure effects is therefore impossible without restrictions on the manner in which $Y_i(\cdot)$ varies with $\bm{d}$. The predominant approach in the literature is to assume interference operates through a low-dimensional vector of sufficient statistics.\footnote{E.g.\ aronow_estimating_2017, basse2019randomization, forastiere2020identification, manski2013identification, toulis2013estimation.} That is, $Y_i(\bm{D})$ is only a function of $\bm{D}$ through a vector-valued {\em exposure mapping}

whose dimension is fixed with respect to $n$, unlike that of $\bm{D}$.

For example, cai2015social run an experiment to study the effect of providing information on the benefits of weather insurance on farmers' take-up of insurance. The authors are interested in spillover effects since a farmer who obtains the information may disseminate it to her social contacts. They estimate linear versions of the model

and $A_{ij}$ is an indicator for whether farmers $i$ and $j$ are friends. Here the exposure mapping $T_i$ is two-dimensional. Variation in its first component identifies the direct effect of the intervention (being offered to attend an information session on weather insurance), while variation in the second identifies a spillover effect.

Exposure mappings, if correctly specified, substantially reduce the dimensionality of the model since we can reparameterize potential outcomes as

where $\mathcal{T}$ is the range of $T(\cdot)$. Interest then centers on “exposure effects”

for $t,t' \in \mathcal{T}$, which measure the average change in potential outcomes in response to counterfactual manipulations of the exposure mapping.

The majority of the literature studies (ref) under the assumption that $T_i$ only depends on treatments assigned to the {\em $K$-neighborhood} of $i$ for some small $K$.\footnote{This is the set of units whose network distance from $i$ is at most $K$, formally defined in (ref).} However, this imposes strong structural restrictions on the underlying outcome process that are incompatible with a variety of models of social interactions studied in the networks literature guilbeault2018complex,jackson2010,manski1993identification. In these models, interference can arise from units outside of the ego's $K$-neighborhood, for any $K$. This is the case for models with endogenous peer effects, where outcomes are functions of the outcomes of neighbors eckles2017design. Another example is the cai2015social setting, where under a simple diffusion model, information obtained by treated units can eventually diffuse to distant alters, which violates (ref).

Important recent work studies {\em misspecified} exposure mappings chin2018central,savje2017average,savje2021causal. The insight of this literature is that standard estimators for (ref) unbiasedly estimate meaningful exposure effects even without imposing (ref) to restrict interference, which indicates a certain robustness of the estimator to more general patterns of interference. However, under what conditions {\em inference} can be made similarly robust is a more challenging question. Without (ref), potential outcomes can depend arbitrarily on the entire assignment vector $\bm{D}$, which makes large-sample inference impossible. These papers accordingly propose a variety of high-level conditions weaker than (ref) that implicitly restrict interference in order to obtain large-sample results. Unfortunately, the connection between these conditions and the literature on social interactions remains unclear. It is an open question whether models in the networks literature violating (ref) (e.g.\ models with endogenous peer effects or other diffusion models) can satisfy these high-level conditions.

\paragraph{Our Contribution.} We study inference on exposure effects with misspecified exposure mappings under a new restriction on interference. We propose a model of {\em approximate neighborhood interference} (ANI), which allows treatments assigned to units further from the ego to have potentially nonzero, but smaller, effects on the ego's response. Unlike the existing literature, we formally verify ANI in well-known models of social interactions. We also show that, under ANI, the data satisfies {\em $\psi$-dependence}, a recently proposed notion of weak network dependence. This enables us to apply limit theorems due to kojevnikov2021limit to establish that, under restrictions on the network topology, standard inverse probability weighting (IPW) estimators are consistent for certain exposure effects and asymptotically normal.

For inference, we consider a network HAC (heteroskedasticity and autocorrelation consistent) variance estimator and characterize its asymptotic bias under a finite population model. We show the bias can be interpreted as the variance of unit-level exposure effects, which is fundamentally unidentified even under no interference. This generalizes the well-known result on conservativeness of the standard variance estimator of the difference-in-means estimate under no interference imbens_causal_2015 to settings with dependence due to interference.

Finally, we propose a novel bandwidth for the network HAC estimator based on the average path length of the observed network. For a given bandwidth rule, the HAC estimator has different rates of convergence depending on whether the average $K$-neighborhood size in $\bm{A}$ grows exponentially or polynomially with $K$. The utility of the average path length is that its magnitude adapts to the neighborhood growth rate to better trade-off bias and variance.

kojevnikov2021limit and kojevnikov2021bootstrap respectively provide consistency results for network HAC and bootstrap variance estimators for $\psi$-dependent network data. Their results pertain to settings in which the data is mean-homogeneous, which often holds in superpopulation models. We extend their results to settings with mean-heterogeneous data, as is the case in finite population models.

choi2017estimation and choi2018using study causal inference without imposing an exposure mapping model. These papers focus on different estimands than ours and assume treatment responses satisfy a monotonicity condition, which we do not require. There is also work on testing for interference, which can be used to test for correct specification of exposure models athey2018exact.

Several papers in econometrics study causal inference under interference baird2018optimal,he2021measuring,lazzati2015treatment,leung2017treatment,vazquez2017identification,viviano2019policy. The second paper studies a dynamic setting, whereas the others focus on a static setting like ours but under correctly specified exposure mappings. Many assume the special case of {\em stratified interference} under which the data consists of many clusters and interference only operates within clusters. We instead consider a single large cluster with known network structure.

The next section states our basic assumptions. In (ref), we present our model of interference (ANI) and large-sample results. We discuss variance estimation in (ref). In (ref), we illustrate the performance of our methods in an empirical application and simulation study. Finally, (ref) concludes. All proofs are given in (ref).

Setup

We consider a finite population model in which the only source of randomness is {\em design uncertainty} abadie2019sampling,imbens2019causal. That is, $\bm{D}$ is the only random quantity. In the special case of no interference, this corresponds to the well-known Neyman causal model. The setup can be viewed as conditioning on the network and potential outcomes, which allows for arbitrary dependence between the two. This allows links to form at higher rates between units with similar unobservables, which corresponds to unobserved homophily, a well-known hindrance to identifying social interactions shalizi2011homophily.

Let $\mathcal{N}_n = \{1, \dots, n\}$ denote the set of units. We assume $\bm{A}$ is an undirected and unweighted network with no self-links, represented as an adjacency matrix with $ij$th entry $A_{ij} \in \{0,1\}$ denoting a potential link between units $i$ and $j$. Let $\mathcal{A}_n$ denote the set of such networks on $n$ units. We assume the components of $\bm{D}$ are independent across units but not necessarily identically distributed, which allows for assignment based on unit covariates and network position. For example, treatment may be assigned “optimally” according to these characteristics viviano2019policy, or randomization may be stratified, as in the empirical application in (ref).

Let $\mathcal{T} \subseteq \mathbb{R}^{d_T}$ be a discrete set. For any $n\in\mathbb{N}$, an {\em exposure mapping} is a function $T \equiv T_n\colon \mathcal{N}_n \times \{0,1\}^n \times \mathcal{A}_n \rightarrow \mathcal{T}$. Most of the literature assumes $T(\cdot)$ is {\em correctly specified} in the sense that (ref) holds, which enables a simple definition of exposure effects (ref). We instead follow the literature on misspecified exposure mappings and employ $T(\cdot)$ only to define useful estimands that summarize treatment and spillover effects but not to restrict the true interference structure. This is a reasonable solution to the task of parsimoniously summarizing the causal effect of a high-dimensional vector $\bm{D}$ on potential outcomes.

For any $n\in\mathbb{N}$ and $i\in\mathcal{N}_n$, a potential outcome $Y_i(\cdot)$ is a mapping from $\{0,1\}^n$ to $\mathbb{R}$. Define the {\em unit-level exposure effect}

and $t,t' \in \mathcal{T}$. This is the difference in unit $i$'s expected response under two different values of the exposure mapping. The estimand of interest is the average effect

also studied by savje2021causal and analogous to estimands proposed by hudgens2008toward but generalized to allow for an incomplete network. We refer to these references for more detailed discussion of interpretation, but the idea is analogous to (ref), which is to compare the average outcomes of units under two different values of the exposure mapping.

Recall $T_i \equiv T(i, \bm{D}, \bm{A})$, and define the generalized propensity score imbens2000role

We estimate $\tau(t,t')$ using the standard IPW estimator, which is unbiased:

Most of the literature focuses on {\em $K$-neighborhood exposure mappings}, requiring $T(\cdot)$ to only be a function of $\bm{d}$ and $\bm{A}$ through $i$'s $K$-{\em neighborhood}, denoted

where $\ell_{\bm{A}}(i,j)$ is the {\em path distance} between $i,j$.\footnote{A {\em path} between $i,j$ is a sequence of links $A_{k_1k_2}, A_{k_2k_3}, \dots, A_{k_{m-1}k_m}=1$ such that $k_1=i$, $k_m=j$, and $k_a \neq k_b$ for all $a,b \in \{1, \dots, m\}$. The {\em length} of this path is $m-1$. The path distance between $i,j$ is the length of the shortest path between them, defined as $\infty$ if $i\neq j$ and no path exists and defined as 0 if $i=j$.} (Note that $i \in \mathcal{N}_{\bm{A}}(i,K)$.) To formalize this requirement, define $\bm{d}_{\mathcal{N}_{\bm{A}}(i,K)} = (d_j\colon j \in \mathcal{N}_{\bm{A}}(i,K))$ and $\bm{A}_{\mathcal{N}_{\bm{A}}(i,K)} = (A_{kl}\colon k,l \in \mathcal{N}_{\bm{A}}(i,K))$, respectively the subvector of $\bm{d}$ and subnetwork of $\bm{A}$ on $\mathcal{N}_{\bm{A}}(i,K)$.

This is a weak restriction on $T(\cdot)$ satisfied by most exposure mappings of interest in the literature. See for example (ref) or the estimands in (ref), where $K=1$.

The next two assumptions impose uniform boundedness on potential outcomes and the generalized propensity score, which will involve bounding these quantities over all $n$. For this to be well-defined, we embed the observed network $\bm{A}$ in a sequence of networks $\{\bm{A}_m\}_{m\in\mathbb{N}}$ where $\bm{A}_m$ is a network on $m$ nodes and $\bm{A}_n = \bm{A}$. We will be imposing conditions on this sequence, which amount to restrictions on the topology of the observed network for large $n$; see (ref) for further discussion.

This requires the generalized propensity score to be uniformly bounded away from 0 and 1 over the sequence of networks. While overlap is standard, it can be restrictive. In particular, it restricts the exposure mapping, network sequence, distribution of treatments, and population. For instance, if $T_i = \mathbf{1}\{\sum_j A_{ij}D_j > 0\}$ and treatments are i.i.d., then overlap holds for $\pi_i(1)$ if treatment assignment is nontrivial and degrees $\sum_j A_{ij}$ are uniformly bounded over $i,n$. However, if $i$'s degree is large, then $\pi_i(1)$ will be close to one since it is highly likely to have at least one of many neighbors treated. At the cost of changing the estimand, overlap can be restored if we instead randomize treatment only to a small subset of “eligible” units and restrict the population to units with an eligible neighbor. Then even with a large degree, the chance of having an {\em eligible} neighbor treated can be far from 0 and 1. This is the design in (ref).

This assumption is standard and can be generalized to uniformly bounded moments.

Approximate Neighborhood Interference

We next present our model of interference. For any $\bm{d} \in \{0,1\}^n$, partition $\bm{d} = (\bm{d}_{\mathcal{N}_{\bm{A}}(i,s)}, \bm{d}_{-\mathcal{N}_{\bm{A}}(i,s)})$, so that $\bm{d}_{-\mathcal{N}_{\bm{A}}(i,s)} = (d_j\colon j \in \mathcal{N}_n\backslash\mathcal{N}_{\bm{A}}(i,s))$. Let $\bm{D}'$ be an independent copy of $\bm{D}$. Define $\bm{D}^{(i,s)} = (\bm{D}_{\mathcal{N}_{\bm{A}}(i,s)}, \bm{D}'_{-\mathcal{N}_{\bm{A}}(i,s)})$, obtained by concatenating the subvector of $\bm{D}$ on $\mathcal{N}_{\bm{A}}(i,s)$ and that of $\bm{D}'$ on $\mathcal{N}_n\backslash\mathcal{N}_{\bm{A}}(i,s)$. Finally, let

This measures interference from “distant” alters, those more than distance $s$ away from the ego. It is the largest expected perturbation of any unit's potential outcome due to redrawing the treatment assignments of distant alters.

This requires interference from distant alters to be negligible for large distances. In the special case of correct specification (ref), there is {\em no} interference from units outside $i$'s $K$-neighborhood, so $Y_i(\bm{D}) - Y_i(\bm{D}^{(i,s)})=0$ for all $s\geq K$ and any $i$. In contrast, ANI allows this difference to be nonzero for all $s$ but requires that it decays with $s$. This means interference from alters beyond a unit's $s$-neighborhood becomes increasingly negligible as we expand the radius $s$. Hence, ANI says that a unit's response is primarily, but not entirely, determined by the assignments of alters close to it.\footnote{(ref) has some similarities with Assumption 6 of chin2018central in bounding the effect of manipulations of treatment assignments of distant units. For a CLT, we do not need a high-level condition analogous to his Assumption 5, which requires correlations between observed outcomes $\{Y_i\}_{i=1}^n$ to be sufficiently weak.}

ANI restricts the network, potential outcomes, and distribution of treatments. The first two are quite evident since the existence of interference is determined by $\bm{A}$ as well as the potential outcomes $Y_i(\cdot)$. For instance, if $\bm{A}$ has no links for any $n$, then ANI holds because for $s>0$, $Y_i(\bm{D}) - Y_i(\bm{D}^{(i,s)})=0$ for any potential outcome model, as no interference is possible without neighbors. On the other hand, for some $Y_i(\cdot)$, ANI can potentially hold for any $\bm{A}$. In the SUTVA case where $Y_i(\bm{d}) = \tilde Y_i(d_i)$ for all $\bm{d}$, this is clearly the case since for any $\bm{A}$ and $s>0$, $Y_i(\bm{D}) - Y_i(\bm{D}^{(i,s)})=0$.\footnote{As pointed out by a referee, since $\theta_{n,s}$ is defined with respect to the design, it may be possible to choose the design to ensure that ANI holds for a given network and potential outcome model. See for example (ref) below, where a sufficient condition for ANI depends on a parameter $\varphi_j$ that is a function of the distribution of $D_j$.}

Social Interactions Models

We next verify ANI for two well-known models of social interactions. Our results yield uniform bounds on $\theta_{n,s}$ that decay exponentially with $s$ under restrictions on the strength of social interactions.

{\bf Linear-in-Means Model.} Consider a network version of the manski1993identification model

where the unobserved heterogeneity $\{\varepsilon_i\}_{i=1}^n$ is uniformly bounded (to impose (ref)) and nonrandom. As usual, to ensure the model is coherent, we assume

The model defines potential outcomes $Y_i(\bm{D})$ through its reduced form

bramoulle2009identification, where $\bm{Y} = (Y_i)_{i=1}^n$, $\bm{\varepsilon}$ is similarly defined, and $\tilde{\bm{A}}$ is the row-normalized version of $\bm{A}$ (divide each row by its sum). The third term roughly says that the impact of treatments assigned to $k$-neighbors is exponentially down-weighted by $\beta^k$. This leads to the following result.

\paragraph{Complex Contagion.} We next consider a model of “complex contagion,” variants of which have been widely studied in the networks literature.\footnote{E.g.\ granovetter1978threshold, guilbeault2018complex, jackson2010, montanari2010spread.} Initialize a dynamic discrete-time process at period 0 at some binary response vector $\bm{Y}^0 \in \{0,1\}^n$, which may be a function of treatments $\bm{D}$ and nonrandom unobserved heterogeneity $\bm{\varepsilon}=(\varepsilon_i)_{i=1}^n$. For some $\mathbb{R}$-valued $\phi(\cdot)$, each unit $i$ at period $t$ updates according to

to obtain new responses $\bm{Y}^t = (Y_i^t)_{i=1}^n$ from last period's responses $\bm{Y}^{t-1}$. The rule says that $i$ selects 1 over 0 if and only if the fraction of neighbors choosing response 1 in the previous period is large enough relative to the heterogeneous threshold $\phi(D_i,\varepsilon_i)$. The parameter $\beta$ measures the strength of social interactions.

Because the setup in (ref) is static, we consider running the dynamic process until the first period $T$ such that $\bm{Y}^T = \bm{Y}^{T-1}$. To ensure such a $T$ exists for any $\bm{Y}^0$, we assume $\beta\geq 0$, the common case of strategic complements milgrom1990rationalizability. We then take $\bm{Y}^T$ as the vector of responses $\bm{Y}$ observed in the data, which yields outcomes $(Y_i(\bm{D}))_{i=1}^n$. Hence, this process implicitly defines potential outcomes.

To verify (ref), we need a condition analogous to (ref), which will be more complicated to state since the model is nonlinear. Define a weighted directed network $\bm{G}$ on $\mathcal{N}_n$ with $ij$th entry $G_{ij} = A_{ij}{\bf E}[\varphi_j]$ for $\varphi_j = \mathbf{1}\{0 < \phi(D_j,\varepsilon_j) \leq \beta\}$. Let

for any $\bar{s} > 0$, so $\lVert\cdot\rVert_\infty$ is the matrix norm induced by the vector $\infty$-norm.

We next discuss the interpretation of $\sup_n \rho_n(\bar{s}) < 1$ in relation to (ref). Let $\bm{I}$ be the $n\times n$ identity matrix. For the linear-in-means model to be coherent, we need $\bm{I}-\beta\tilde{\bm{A}}$ to be invertible, which is true provided $\lvert\beta\rvert<1$. Since $\tilde{\bm{A}}$ is row-normalized, this is equivalent to

where $\lambda_\text{max}(\cdot)$ is the spectral radius bramoulle2009identification. On the other hand, $\lVert\bm{G}^s\rVert_\infty^{1/s} \stackrel{s\rightarrow\infty}\longrightarrow \lambda_\text{max}(\bm{G})$ by Gelfand's formula, so for any $\epsilon>0$, we can choose $\bar{s}$ large enough such that

xl2015. The left-hand side is clearly analogous to (ref). The difference is that, in $\bm{G}$, we weight each potential link $A_{ij}$ by ${\bf E}[\varphi_j]$, whereas in $\beta\tilde{\bm{A}}$, the weight is $\beta/\sum_k A_{ik}$. Both weights are monotonically increasing in $\beta$, so both (ref) and (ref) restrict the strength of social interactions.

Weak Dependence

Define $Z_i = (\bm{1}_i(t) \pi_i(t)^{-1} - \bm{1}_i(t') \pi_i(t')^{-1}) Y_i$, so that $\hat\tau(t,t') = n^{-1}\sum_{i=1}^n Z_i$. For large-sample inference, we would like the data $\{Z_i\}_{i=1}^n$ to be at most weakly dependent. Recall that $\ell_{\bm{A}}(i,j)$ is the path distance between $i,j$ in $\bm{A}$. Since treatments are independent, the indicators are weakly dependent in the sense that $\bm{1}_i(t) \perp\!\!\!\perp \bm{1}_j(t)$ if $\ell_{\bm{A}}(i,j) > 2K$ by (ref). At first glance, $\{Y_i\}_{i=1}^n$ seems to be strongly dependent since $Y_i$ is a function of $\bm{D}$ for all $i$. However, ANI requires $Y_i$ to primarily depend on the treatments of nearby units, which suggests that distant units have weakly dependent outcomes. We next formalize this idea.

We first define a notion of weak network dependence due to kojevnikov2021limit. For any $H,H' \subseteq \mathcal{N}_n$, define $\ell_{\bm{A}}(H,H') = \min\{\ell_{\bm{A}}(i,j)\colon i \in H, j \in H'\}$. Let $\bm{Z}_H = (Z_i\colon i \in H)$, $\mathcal{L}_d$ be the set of bounded, $\mathbb{R}$-valued, Lipschitz functions on $\mathbb{R}^d$, and

This extends temporal $\psi$-dependence doukhan1999new to network data by using path in place of temporal distance. The concept says two sets of observations $\bm{Z}_H$ and $\bm{Z}_{H'}$ have small covariance if they are sufficiently distant.

Let $i^*(\cdot)$ be the identity function $x \mapsto x$ on $\mathbb{R}$ and $\mathcal{Z} \subseteq \mathbb{R}$ be any compact set such that $Z_i \in \mathcal{Z}$ for any $n\in\mathbb{N}$ and $i\in\mathcal{N}_n$. Such a set exists by Assumptions (ref) and (ref). For any Lipschitz function $f\colon\mathbb{R}^d\rightarrow\mathbb{R}$, let $\text{Lip}(f)$ its Lipschitz constant and $\lVertf\rVert_\infty=\sup_{x\in\mathcal{Z}^d} \lvertf(x)\rvert$. Let $K$ be the constant in (ref) and $\lfloor s \rfloor$ be $s$ rounded down to the nearest integer.

Large-Sample Theory

In what follows, define $\tilde\theta_{n,s}$ as in (ref). Having established that $\{Z_i\}_{i=1}^n$ is $\psi$-dependent, we can apply results due to kojevnikov2021limit to show that $\hat\tau(t,t')$ is consistent and asymptotically normal. Their results require $\tilde\theta_{n,s}$ to decay to zero fast enough, and the speed of decay depends on the network topology, in particular the growth rate of $s$-neighborhood sizes. Intuitively, ANI says that $Z_i$ depends primarily on units in $\mathcal{N}_{\bm{A}}(i,s)$. Hence, if the typical size of these neighborhoods grows rapidly with $s$, then weak dependence requires that this be counterbalanced by having the covariances ($\tilde\theta_{n,s}$) decay to zero faster with $s$.

The next assumptions formalize these ideas. To clarify their interpretation, in (ref), we verify them for networks with polynomial and exponential neighborhood growth rates when $\tilde\theta_{n,s}$ decays exponentially with $s$, as in the examples in (ref).

{\bf Limit Sequence.} Our results take $n\rightarrow\infty$ along a sequence of networks $\{\bm{A}_n\}_{n\in\mathbb{N}}$ defined prior to (ref). The results hold for any sequence satisfying our assumptions. The design may implicitly depend on $n$, so the distribution of treatments may also change along the sequence, so long as treatments remain independent across units. Additionally, potential outcomes $Y_i(\cdot)$ and the exposure mapping $T(\cdot)$ may vary with $n$ (and must do so due to the dimensions of their arguments), so long as they satisfy Assumptions (ref) and (ref). We emphasize that our results assume $K$ and $\mathcal{T}$ do {\em not} depend on $n$, but all other quantities may vary along the sequence, within the confines of the stated assumptions.

Let $\mathcal{N}^\partial_{\bm{A}}(i,s) = \{j\in\mathcal{N}_n\colon \ell_{\bm{A}}(i,j)=s\}$ be the {\em $s$-neighborhood boundary} of $i$, the set of units exactly distance $s$ from $i$, and $M_n^\partial(s) = n^{-1} \sum_{i=1}^n \lvert\mathcal{N}^\partial_{\bm{A}}(i,s)\rvert$, its average size.

This corresponds to Assumption 3.2 of kojevnikov2021limit. It restricts the network topology through $M_n^\partial(s)$ and the degree of interference through $\tilde\theta_{n,s}$. Since the former grows with $s$, the latter must decay to zero faster for the sum to be $o(n)$. Furthermore, since $\tilde\theta_{n,s}=1$ for $s\leq 2\max\{K,1\}$, this implies $M_n^\partial(1) = o(n)$, which is a restriction on network density that rules out, for instance, $\bm{A}$ being complete.

It is useful to compare this to its analog for $\alpha$-mixing spatial processes. Consider, for example, Assumption 3(b) of jenish2009central, which essentially requires $\sum_{s=1}^\infty s^{d-1} \alpha(s) < \infty$, where $d$ is the dimension of the underlying space and $\alpha(s)$ is the $\alpha$-mixing coefficient, which measures dependence between sets of observations at spatial distance $s$ apart. In the spatial setting, the $s$-neighborhood boundary of $i$ is the set of units at any distance $h \in [s,s+1)$ from $i$. By their Lemma A.1(iii), the size of this set is $O(s^{d-1})$. Thus, we have an analogous trade-off between the sizes of spatial $s$-neighborhood boundaries and the rate of decay of the mixing coefficient.

Inspection of the proof shows that we can sharpen the result to $\lvert\hat\tau(t,t') - \tau(t,t')\rvert = O_p(n^{-1/2})$ if we strengthen (ref) to $\sum_{s=0}^n M_n^\partial(s) \tilde\theta_{n,s} = O(1)$.

Asymptotic normality requires a stronger version of (ref). Let $M_n(s, k) = n^{-1}\sum_{i=1}^n \lvert\mathcal{N}_{\bm{A}}(i, s)\rvert^k$, the $k$th moment of the $s$-neighborhood size, and

This is the set of paired couples $(i,j)$ and $(k,l)$ such that the units within each couple are at most path distance $m$ apart from one another, and the two pairs are exactly path distance $s$ apart. Define $\sigma_n^2 = \text{Var}(n^{-1/2}\sum_{i=1}^n Z_i)$.

This strengthens (ref) and corresponds to Assumption 3.4 of kojevnikov2019limit.\footnote{The publication version of their paper kojevnikov2021limit formulates this assumption slightly differently.} Similar to (ref), the first term in (ref) requires $\tilde\theta_{n,s}$ to decay to zero fast enough relative to $s$-neighborhood sizes. The second term restricts $s$-neighborhood growth rates, while the third requires sufficiently fast decay of $\tilde\theta_{n,s}$. See (ref) for a discussion of the plausibility of this assumption.

Variance Estimation

For large-sample inference, we consider the variance estimator

where $Z_i$ is defined at the top of (ref) and $b_n\geq 0$ is a bandwidth parameter discussed below in (ref). When $b_n=0$, this reduces to the sample variance of the $Z_i$'s, which is only a valid estimator under no interference. Choosing $b_n>0$ places nonzero weight on pairs at most $b_n$ apart in the network, which accounts for possible autocorrelation.\footnote{The estimator is simple to compute. First calculate the path distance matrix $(\ell_{\bm{A}}(i,j))_{i,j\in\mathcal{N}_n}$, which can be done very efficiently for sparse networks using Dijkstra's algorithm (e.g.\ the Python function {\tt dijkstra} in the {\tt scipy.sparse.csgraph} module). Then for the matrix $\bm{P} = (\mathbf{1}\{\ell_{\bm{A}}(i,j)\leq b_n\})_{i,j\in\mathcal{N}_n}$ and vector $\tilde{\bm{Z}} = (n^{-1/2} (Z_i - \hat\tau(t,t')))_{i=1}^n$, we have $\hat\sigma^2 = \tilde{\bm{Z}}'\bm{P}\tilde{\bm{Z}}$.}

Under correctly specified exposure mappings, we can choose $b_n = 2K$, as in leung2017treatment, since $\bm{1}_i(t) \perp\!\!\!\perp \bm{1}_j(t)$ if $\ell_{\bm{A}}(i,j) > 2K$. For misspecified exposure mappings, $b_n$ must grow with $n$ at a rate depending on the network topology. In this case, (ref) corresponds to a HAC estimator familiar from the time series literature but using network distance in place of temporal distance. This estimator has been previously used in practice acemoglu2015state, and its formal properties were first studied by kojevnikov2021limit in a superpopulation setting. Building on their work, we characterize its behavior in a finite population model.

\paragraph{Choice of Bandwidth.} For $\hat\sigma^2$ to have good large-sample properties, we need to restrict the rate at which $\lvert\mathcal{N}_{\bm{A}}(i,b_n)\rvert$ diverges with $n$ (see (ref)). That is, how fast $b_n$ can diverge depends on how rapidly $K$-neighborhood sizes grow with $K$. In spatial settings, the rate of growth is polynomial in $K$, so $b_n$ is allowed to diverge at a polynomial rate. (A faster rate is better for bias but worse for variance.) However, in network settings, the rate can be exponential.

Based on the analysis in (ref), we suggest $b_n$ be chosen as follows. Let $\delta(\bm{A}) = n^{-1} \sum_{i,j} A_{ij}$ be the average degree, $\mathcal{L}(\bm{A})$ the {\em average path length} (APL), and

where $\lfloor\cdot\rceil$ means round to the nearest integer.\footnote{The formula assumes $\delta(\bm{A})>1$, which is typical in practice. The APL is the average value of $\ell_{\bm{A}}(i,j)$ over all pairs in the largest component of $\bm{A}$. A {\em component} of a network is a connected subnetwork such that all units in the subnetwork are disconnected from those not in the subnetwork.} To account for correlation in $\{\bm{1}_i(t)\}_{i=1}^n$ discussed in (ref), we set $b_n$ at least equal to $2K$. The fractions $1/2$ and $1/3$ are due to (ref) below (see (ref)). We suggest in practice the researcher report results for several bandwidths in a neighborhood of (ref).

The purpose of comparing $\mathcal{L}(\bm{A})$ and $\log n / \log \delta(\bm{A})$ is to determine whether $K$-neighborhood sizes grow approximately exponentially or polynomially with $K$. As discussed in (ref), in the exponential case, the difference between these two statistics typically converges to zero, whereas in the polynomial case, $\mathcal{L}(\bm{A})$ is much larger, having polynomial order. See, for example, the simulations in (ref), which show at least a four-fold difference in APL between the two regimes. Thus, (ref) selects a bandwidth of logarithmic (polynomial) order when neighborhood growth rates are approximately exponential (polynomial).\footnote{In the exponential case, we need $b_n = O(\log n)$ for $\text{Var}(\hat\sigma^2)$ to be small, which (ref) accomplishes since $\mathcal{L}(\bm{A}) \approx \log n / \log \delta(\bm{A})$ in this regime. In the polynomial case, $b_n = O(\log n)$ is also valid, but the bias then vanishes at an extremely slow rate kojevnikov2021limit. Our choice of $\mathcal{L}(\bm{A})^{1/3}$ substantially improves this rate since $\mathcal{L}(\bm{A})$ is then polynomial in $n$.}

\paragraph{Bias of $\hat\sigma^2$.} Define

The former is an “oracle” version of $\hat\sigma^2$ that replaces $\hat\tau(t,t')$ with $\tau_i(t,t')$, while the latter is a bias term. (ref) below establishes that

Equation (ref) says that the oracle estimator is consistent for the variance, and (ref) says that our estimator is biased. The source of bias is mean-heterogeneity: $\hat\tau(t,t')$ is consistent for $\tau(t,t')$ but not $\tau_i(t,t')$, which is heterogeneous across units.

The bias $R_n$ has the form of a HAC estimate of the variance of the unit-level exposure effects. It is helpful to compare this to the case of no interference, where $T_i = D_i$, $t=1$, and $t'=0$, so that $\tau(t,t')$ is the usual average treatment effect (ATE). Knowing that units are independent, we can choose $b_n=0$, in which case

This is the well-known asymptotic bias of the standard variance estimator for the difference-in-means estimate of the ATE imbens_causal_2015. It measures the variance of the unit-level treatment effects and is generally impossible to estimate in the finite population setting, so the variance estimator is conservative. In the special case of homogeneous unit-level treatment effects, meaning $\tau_i(t,t')$ does not vary with $i$, the bias is zero, a property also shared by our $R_n$. Thus, (ref) generalizes Neyman's well-known result on conservative variance estimation to a setting with interference. The additional covariance terms in $R_n$ weighted by $\mathbf{1}\{\ell_{\bm{A}}(i,j)\leq b_n\}$ account for dependence due to interference.

To show consistency of $\hat\sigma^2$, define

This is similar to, and evidently contains, $\mathcal{H}_n(s,m)$ from (ref).

In (ref), we use these conditions to derive (ref) and illustrate their plausibility for some general classes of graphs. Part (a) strengthens (ref), while (b)--(d) regulate $b_n$. Part (b) allows us to replace $\hat\tau(t,t')$ in $\hat\sigma^2$ with its expectation. Part (d) is used to derive the asymptotic bias. It is very similar to the first requirement of (ref), with $b_n$ and $\mathcal{J}_n(s,\cdot)$ in place of $m_n$ and $\mathcal{H}_n(s,\cdot)$, respectively.

Numerical Illustrations

Empirical Application

We revisit a network experiment analyzed in paluck2016changing and aronow_estimating_2017 that studies the effect of an anti-conflict intervention on adolescent social norms for antagonistic behavior, including harassment, rumor-mongering, social exclusion, and bullying. In the experimental design, 28 of 56 schools are first randomized into treatment. Then within treated schools, a subset of students are selected as eligible for treatment based on covariates, and half of eligibles are block-randomized into treatment. Treated students are invited to participate in bi-monthly meetings that follow an anti-conflict curriculum designed in part by the researchers of the study. At these meetings, a trained adult leader helps students identify social conflicts at their school and design strategies to reduce conflict.

aronow_estimating_2017 and part of the analysis of paluck2016changing examine the causal effect of the offer to participate on endorsement of anti-conflict norms. This is measured by self-reports of wearing a wristband disseminated as part of the program as a reward to students observed engaging in conflict-mitigating behavior. Through the course of the experiment, over 2500 wristbands were disseminated and tracked. We study exposure effects similar to those used in the application of aronow_estimating_2017. Unlike their analysis, we restrict the data to the five largest treated schools to illustrate what can be learned from data on a few large networks. In each of our schools, the number of eligibles is exactly 64.

We estimate a treatment and a spillover effect. For the latter, the exposure mapping is $T_i = \bm{1}\{\sum_j A_{ij}D_j > 0\}$, an indicator for whether at least one friend is offered treatment. As in aronow_estimating_2017, to ensure overlap, we restrict to the “spillover population” consisting of students that have at least one eligible friend ($n=1685$). For the treatment effect, the exposure mapping is $T_i = D_i$, and we restrict to the “treatment population” consisting of students eligible for treatment ($n=320$). For both, we compute path distances and average degree in the variance estimator (discussed below) using the full network.

Networks are measured by asking students to name up to ten students at the school “whom they chose to spend time with in the last few weeks, either in school, out of school, or online.” Consequently, $\bm{A}$ is directed. When computing the number of treated friends for the exposure mappings, we use the directionality of links. However, when computing network neighborhoods for our variance estimator, we ignore the directionality of links to conservatively define larger neighborhoods and avoid taking a stance on neighborhood definitions for directed networks.

We next provide some summary statistics. Within the treatment population, the average outcome $Y_i$ is 0.16 (SD 0.37), and by block randomization, exactly 50 percent are treated. Within the spillover population, the average outcome is 0.11 (SD 0.32), and 58 percent (SD 0.49) have at least one treated friend. The data includes the blocks in which eligible students are block-randomized, so we can compute the propensity scores $\pi_i(t)$ for each student using the hypergeometric distribution. For the exposure mapping $T_i = \bm{1}\{\sum_j A_{ij}D_j > 0\}$, given $N$ eligible neighbors, $\pi_i(0)$ is the chance of having 0 out of $N$ successes when drawing without replacement. We find $n^{-1} \sum_{i=1}^n \pi_i(1) = 0.597$, which is close to empirical proportion of 58 percent.

The average out-degree $n^{-1} \sum_{i,j} A_{ij}$ is 7.96. The APL is small, on average 3.37 across our five schools. Since there are $n=3306$ students, $\log n / \log \delta(\bm{A}) = 3.96$, which is very close to 3.37. Thus, given $K=1$, our suggested bandwidth (ref) is $b_n=2$. We report results for the range of bandwidths $\{0,\dots,3\}$, noting that $0$ corresponds to standard errors in the no-interference case.

(ref) presents the results. The first row is the IPW estimator for the indicated exposure effect, and the last four rows are standard errors for the indicated bandwidths. We find a large treatment effect of 0.15, which is significant at the 5 percent level across all bandwidths. The spillover effect is smaller at 0.04, with larger standard errors, and is statistically insignificant for our suggested bandwidth $b_n=2$ at the 5 percent level. The small spillover estimate is largely in line with the estimates implied by Figure 3C of paluck2016changing. It does not contradict the overall message of their paper since they find, for example, sizeable spillover effects when comparing treated and untreated schools. In contrast, our analysis so far only makes comparisons within treated schools, using only a subsample of five schools.

To compare treated and untreated schools, note that for the latter, $Y_i=0$ for all $i$ by design. Then our estimate of $\hat\mu(0)$ for the treatment effect shows that even untreated units are 8.8 percentage points more likely to wear wristbands in treated compared to untreated schools. We obtain standard errors for $\hat\mu(0)$ by replacing $Z_i$ and $\hat\tau(t,t')$ in (ref) with $\bm{1}_i(t) \pi_i(t)^{-1} Y_i$ and $\hat\mu(0)$. For $b_n = 0,\dots,3$, the standard errors range from $0.012$ to $0.017$, all of which imply a statistically significant effect. Overall, these results indicate that, despite the potential conservativeness of our estimator due to the bias term $R_n$, they can still deliver reasonable standard errors.

Monte Carlo

To study the finite sample properties of our estimators, we simulate data from the two response models studied in (ref) using two models of network formation calibrated to the school data from (ref). For the linear-in-means model, $Y_i = V_i(\bm{D},\bm{A},\bm{\varepsilon})$ for

and $(\alpha,\beta,\delta,\gamma) = (-1,0.8,1,1)$. For the complex contagion model, we set $Y_i = \bm{1}\{V_i(\bm{D},\bm{A},\bm{\varepsilon}) > 0\}$ and $(\alpha,\beta,\delta,\gamma) = (-1,1.5,1,1)$.

We simulate $\bm{A}$ from configuration and random geometric graph (RGG) models. The former is calibrated to the empirical out-degree sequence $(\sum_{j=1}^n A_{ij})_{i=1}^n$ of the schools used in (ref). This model (approximately) draws an undirected network uniformly at random from the set of all networks with this degree sequence jackson2010. An RGG is a spatial network where units only link with geographically close alters: $A_{ij} = \bm{1}\{\lVert\rho_i-\rho_j\rVert \leq r_n\}$ for $\rho_i \stackrel{iid} \sim \mathcal{U}([0,1]^2)$ and $r_n = (\kappa/(\pi n))^2$. Since $\kappa$ is the limiting expected degree of the model penrose2003, we set it equal to the average of the empirical out-degree sequence of the schools.

Let $\{\nu_i\}_{i=1}^n \stackrel{iid}\sim \mathcal{N}(0,1)$ be independent of $\bm{A}$. For the configuration model, we take $\varepsilon_i = \nu_i$. For the RGG, we instead use $\varepsilon_i = (\rho_{i1}-0.5) + \nu_i$, adding the centered first component of $i$'s “location” $\rho_{i1}$. Since units with similar $\rho_{i1}$'s are more likely to form links, this generates unobserved homophily.

Both models yield networks with approximately the same average degree of about 8. For our purposes, the main distinction between them is the APL. The configuration model is theoretically known to have a small APL of logarithmic order in the network size van2016random, while the RGG model seems to generate a much larger APL of polynomial order friedrich2013diameter. We choose different bandwidths for the two according to (ref).

Note that the value of the peer effect $\beta$ in both outcome models is actually quite large. The calculations in Appendices (ref) and (ref) suggest that it is larger than what our weak dependence conditions actually require under the configuration model (but not the RGG) since it generates a low APL (and hence approximately exponential neighborhood growth rates). Of course, these conditions are likely stronger than necessary, which may explain the good performance in both designs below.

To show different population sizes, we report results using the largest, two largest, and four largest of the treated schools when calibrating the network models. In all cases, we pool the degree sequences across the schools to treat them as one single network. Following the design in (ref), we randomly assign treatments to only units classified as eligible in the data with probability 0.5.

We compute estimates and standard errors for the spillover effect $\tau(1,0)$ in (ref), whose exposure mapping is $T_i = \bm{1}\{\sum_j A_{ij}D_j > 0\}$, again restricting to the population of units with an eligible neighbor. The standard errors use the bandwidth (ref). Given the APL differences discussed above, which can also be seen in Tables (ref) and (ref), the RGG uses $\tilde b_n = \mathcal{L}(\bm{A})^{1/3}$, while the configuration model uses $\tilde b_n = \mathcal{L}(\bm{A})/2$.

To illustrate the performance of our bandwidth rule across different networks, we redraw $\bm{D}$, $\bm{A}$, and $\bm{\varepsilon}$ for each of the 10k simulation draws. This corresponds to a superpopulation design, so we expect our standard errors to yield confidence intervals with coverage close to the nominal rate of 0.95. We report “oracle” standard errors, which correspond to $\text{Var}(\hat\tau(1,0))^{1/2}$, approximated by taking the standard deviation of $\hat\tau(1,0)$ over 10k separate simulation draws. We also report “naive” i.i.d.\ standard errors to illustrate the degree of dependence in the data.

Tables (ref) and (ref) present results for the configuration and RGG models, respectively. “Our SE” displays the standard error obtained from our variance estimator and “Our Coverage” the empirical coverage of the corresponding two-sided 95% confidence interval. The effective sample size of $\hat\mu(t)$ is $\hat n(t) = \sum_i \bm{1}_i(t)$. “Network Size” is the number of units in the network, whereas the sample size is $\hat n(1) + \hat n(0)$ since, as in the empirical application, we only use units with an eligible neighbor for overlap.

The results show that the oracle coverage rates are all very close to 0.95, which illustrates the quality of the normal approximation. Our standard errors perform well, with coverage quite close to 0.95. This is in spite of fairly large peer effects $\beta$, which can be seen in the magnitudes of the spillover effect estimates. By contrast, naive standard errors are severely anti-conservative.

Conclusion

The literature on causal inference under interference typically assumes $K$-neighborhood exposure mappings are correctly specified. This implies a limited model of interference in which units further than distance $K$ from the ego have no effect on the ego's response for some small, known $K$. Such a model is incompatible with those studied in the large theoretical and empirical literature on social interactions, in which units arbitrarily far from the ego can influence the ego's outcome. This paper proposes a richer model of “approximate neighborhood interference” (ANI) under which treatments assigned to distant alters have a nonzero effect on the ego's outcome, but the effect declines with network distance. Unlike with existing models of interference, we show that ANI is satisfied by well-known models of social interactions. We also show that ANI is useful for large-sample inference, proving that IPW estimators are consistent and asymptotically normal.

We consider a network HAC variance estimator for inference robust to misspecification of exposure mappings. We show the estimator is biased in a finite population setting. The bias term captures the variance of the unit-level treatment effects, which generalizes the well-known result on conservative variance estimation under no interference to settings with dependence due to interference. Finally, we propose a new bandwidth for the HAC estimator, which trades off bias and variance in a manner that adapts to the growth rate of $K$-neighborhoods.