Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
91,601 characters · 18 sections · 143 citation commands
Causal inference in network experiments: regression-based analysis and design-based properties
\onehalfspacing
\thispagestyle{empty}
Keywords: Covariate adjustment, exposure mapping, interference, model misspecification, network-robust standard error, weighted least squares.
Network experiments have gained growing interest across various fields, including economics, social science, public health, and tech companies Jackson2008, Valente2010, BlakeCoey2014, AngelucciDiMaro2016, Aral2016, Breza2016, AtheyImbens2017, AtheyEcklesImbens2018a, AronowEcklesSamii2021. They present an exceptional avenue to delve into the intricacies of interactions among units. Important examples of such experiments include Sacerdote2001, MiguelKremer2004, BandieraRasul2006, BakshyRosennMarlow2012, BanerjeeChandrasekharDuflo2013, BursztynEdererFerman2014, CaiJanvrySadoulet2015, PaluckShepherdAronow2016, BeamanDillon2018, HaushoferShapiro2018, and CarterLaajajYang2021. These experiments transcend the conventional framework of individual-level randomization by exploring the effects of treatments not only on the treated individuals but also on their peers. This introduces the concept of “interference,” which challenges the “stable unit treatment value assumption” (SUTVA) that rules out interference in classic causal inference.
Over the last decade, the study of social interactions and peer effects through structural models has gained considerable attention Manski1993, Graham2008, BramoulleDjebbariFortin2009a, Goldsmith-PinkhamImbens2013. Distinguishing between the influence of peers' outcomes (endogenous peer effects) and the influence of peers' characteristics (contextual peer effects) can become challenging due to the simultaneous behavior of interacting agents. This challenge is known as the “reflection problem” Manski1993. Angrist2014 criticized several econometric approaches to estimating peer effects. Without covariates, outcome-on-outcome regressions either reflect a tautological identity or capture group-level clustering without behavioral meaning. With covariates, the resulting estimates may be biased by measurement error and other factors, leading to spurious evidence of peer effects. Paula2017 and BramoulleDjebbariFortin2020 relate the counterexample proposed in Angrist2014 to a well-known instance of non-generic identification failure, initially noted by Manski1993 and also demonstrated by BramoulleDjebbariFortin2009a.
An expanding volume of literature explores scenarios with interference of arbitrary but known forms which in turn requires researchers to make specific assumptions about the extent of interference. Many papers assume correctly specified exposure mappings for inference AronowSamii2017, BairdBohrenMcIntosh2018, Vazquez-Bare2022, Owusu2023. These mappings impose assumptions on the interference structure in the experiment, where the treatment assignment vector affects potential outcomes through a low-dimensional function Manski2013, AronowSamii2017. This approach can be critiqued for typically ruling out endogenous peer effects. Some other papers assume “partial interference” Sobel2006, HudgensHalloran2008, UganderKarrerBackstrom2013, KangImbens2016, LiuHudgensBecker-Dreps2016, BasseFeller2018, QuXiongLiu2021, AlzubaidiHiggins2023, where units are partitioned into separate clusters, and interference is restricted to occur exclusively among units within the same cluster. Conversely, more recent literature further relaxes the partial interference assumption and studies interference of general forms SAVJEARONOWHUDGENS2021, Viviano2023.
Leung2022 proposed to estimate exposure effects under “approximate neighborhood interference” (ANI) while allowing for misspecification of exposure mappings. ANI refers to the situation where treatments assigned to individuals further from the focal unit have a smaller, but potentially nonzero, effect on the focal unit's response. Leung2022 verified that ANI is applicable to well-known models of social interactions, such as the network version of the linear-in-means model Manski1993 and the complex contagion model Granovetter1978, both of which allow for endogenous peer effects. He considered the Horvitz--Thompson estimator and studied its consistency and asymptotic normality. For inference, he proposed a network Heteroskedasticity and Autocorrelation Consistent (HAC) covariance estimator, and studied its asymptotic bias for estimating the true covariance. However, he did not derive the point and covariance estimator directly from regression-based analysis, which is our focus.
Our paper builds upon Leung2022, which accommodates a single large network. We enrich the discussion of the regression estimators from the design-based perspective, with a special emphasis on network experiments. The design-based inference makes weak distributional assumptions about outcome models and relies solely on the randomization mechanism. We focus on the Hájek estimator, which is numerically identical to the coefficient derived from the weighted-least-squares (WLS) fit involving unit data that relies on the inverse probability of exposure mappings AronowSamii2017. The regression-based approach offers three notable advantages. First, it is easy to implement without too much additional programming. Second, it can provide standard errors through the same WLS fit. Third, it allows for incorporating covariates into the analysis, which can potentially increase the estimation precision if the covariates are predictive of the outcome. Moreover, we examine the asymptotic performance of the regression-based network HAC estimator and prove results that justify the regression-based inference for network experiments from the design-based perspective.
Unlike their spatial or time-series counterparts, network HAC estimators lack a theoretical guarantee of positive semi-definiteness Kojevnikov2021. Moreover, they are known to have poor finite-sample properties Matyas1999. In network experiments, the asymptotic bias of the HAC estimator can be negative under interference, resulting in undercoverage of the associated confidence interval. To address these concerns, we propose a modified HAC estimator that ensures positive semi-definiteness and asymptotic conservativeness, which also performs well in finite-sample simulation.
Furthermore, we delve into the subject of covariate adjustment. Proper covariate adjustment can enhance the accuracy of estimators in randomized experiments by accounting for the imbalance in pretreatment covariates. Recall the results in the classical completely randomized treatment-control experiment. The regression framework offers a versatile approach to incorporating covariate information with a potential of enhancing asymptotic efficiency by including the interactions of the treatment and covariates Fisher1935, Lin2013, NegiWooldridge2021. An expanding body of literature explores the design-based justification of regression-based covariate adjustment with different types of experimental data Fogarty2018, ChangMiddletonAronow2021, SuDing2021, ZhaoDing2022, WangSusukidaMojtabai2023, ZhaoDingLi2024. Our paper studies the theoretical properties of covariate adjustment in network experiments and demonstrates the potential efficiency gain in simulation and empirical application under reasonable data-generating processes.
\paragraph*{Organization of the paper} Section (ref) sets up the framework for the design-based inference in network experiments, reviews the Horvitz--Thompson and Hájek estimators, and introduces the main assumptions from Leung2022. Section (ref) reviews the Hájek estimator recovered from the WLS fit AronowSamii2017, proposes the regression-based HAC covariance estimator, and analyzes its asymptotic bias. Because the covariance estimator can be anti-conservative, we propose a modified, positive semi-definite covariance estimator. Section (ref) considers additive and fully-interacted covariate adjustment to the WLS fit, describes associated asymptotic properties, proposes modified covariance estimators, and studies their asymptotic properties. Section (ref) studies the finite-sample performance of our point and covariance estimators based on simulation and illustrates the practical relevance of our results by re-analyzing the network experiment in PaluckShepherdAronow2016. Section (ref) discusses the extension to continuous exposure mappings. The appendix includes all the proofs and intermediate results.
\paragraph*{Notation} Let $\mathbb{N}$ denote the set of all non-negative integers. Let ${I}_{m}$ be an $m \times m$ identity matrix and $\iota_{m}$ be an $m \times 1$ vector of ones. We suppress the dimension $m$ when it is clear from the context. Unless stated otherwise, all vectors are column vectors. Let $1(\cdot)$ be the indicator function. Let $\|\cdot\|$ denote the Euclidean norm, i.e., $\|w\| = \sqrt{w^\top w}$ for $w\in\mathbb{R}^v$. The terms “regression” and “HAC covariance” refer to the numerical outputs of the WLS fit without any modeling assumptions; we evaluate their properties under the design-based framework. We use “IID” and “CLT” to denote “independent and identically distributed” and “central limit theorem,” respectively.
We consider a finite population model that conditions on potential outcomes and networks while viewing the treatment assignment as the sole source of randomness. This approach follows the design-based framework ImbensRubin2015, AronowSamii2017, LiDing2020a, AbadieAtheyImbens2020, Leung2022, Chang2023. Let $\mathcal{N}_n = \{1,\ldots,n\}$ denote the set of units. The network structure is undirected, unweighted, has no self-links, and can be described using an adjacency matrix $ {A}=(A_{ij})_{i,j=1}^n$ with the $(i, j)$th entry $A_{ij} \in \{0, 1\}$ indicating the connection between units $i$ and $j$. Let $\mathcal{A}_n$ denote the set of all possible networks with $n$ units. The assignment of treatments is represented by a binary vector $ {D} = (D_i)^n_{i=1}$, where each $D_i$ is a binary variable indicating whether unit $i$ has been assigned to the treatment.
We define the potential outcome for each unit $i$ as $Y_i( {d})$, which represents the outcome of unit $i$ under the hypothetical scenario in which the units on the entire network are assigned the treatment vector ${d} = (d_i)^n_{i=1}\in \{0,1\}^n$. From the notation, $Y_i( {d})$ depends not only on $d_i$, the treatment assignment of unit $i$, but also on the treatment assignments of all other units. This results in “interference” or “spillover” between units, which is not accounted for in the standard potential outcomes model under SUTVA. We adopt the design-based framework in which the potential outcomes $Y_i(d)$'s and network $A$ are fixed, whereas the distribution of $D$ is known and does not depend on $Y_i(d)$'s and $A$.
With binary treatment $D_i$'s, we have $2^n$ potential outcomes for each unit. We utilize the exposure mapping as defined by AronowSamii2017 or the effective treatment mapping introduced by Manski2013 for dimensionality reduction. For any $n$, an exposure mapping is a function $T: \mathcal{N}_n \times\{0,1\}^n \times \mathcal{A}_n \rightarrow \mathcal{T}$, which maps the units, the treatment assignment vector, and the network structure to exposures received by a unit. We focus on the regime in which $|\mathcal{T}|$ is finite and fixed. With correctly specified exposure mapping, we can simplify the potential outcomes as $Y_i(d) = {Y}_i(t)$ because $Y_i(D)$ depends on $D$ only through $T_i = T(i, D, A)$. We follow Leung2022's framework and allow for misspecified exposure mappings: $T_i$ need not correctly capture how others affect an individual's potential outcome. With the potential outcomes $Y_i(d)$, we can define the unit $i$'s expected response under exposure mapping value $t$ as
which equals the expected potential outcome of unit $i$ over all possible treatment assignment vectors given the exposure mapping value at $t$. Let $\mu(t) = n^{-1} \sum_{i=1}^n \mu_i(t)$ be the finite-population average and $ {\mu} = (\mu(t):t\in\mathcal{T})$ be the $|\mathcal{T}|\times 1$ vector containing all the $\mu(t)$'s corresponding to exposure mapping values $t\in\mathcal{T}$. We will discuss inference of the general estimand $\tau = G {\mu}$, where $G$ is an arbitrary contrast matrix, and the key lies in estimating $ {\mu}$. We focus on estimators of the form $\hat{\tau} = G \hat{Y}$, where $\hat{Y}$ is some regression estimator of $\mu$. Although we focus on regression-based point estimators and standard errors, our theory holds under the design-based framework, which assumes that the randomness comes solely from the design of network experiments and allows for misspecification of the regression models.
While the theory can accommodate misspecified exposure mappings, this flexibility comes at the cost of complicating the causal interpretation. When the exposure mapping is correctly specified, we have $Y_i(d) = {Y}_i(t)$ for $t \in \mathcal{T}$, allowing the average expected response to simplify to $\mu(t) = n^{-1}\sum_{i=1}^n {Y}_i(t)$. In this case, the estimand ${\tau} = G \mu$ becomes independent of the treatment assignment and has a clear causal interpretation. However, when the exposure mapping is misspecified, $\mu(t)$ represents a weighted average of all potential outcomes, where the weights correspond to $\mathbb{P}\left({D}={d} \mid T_i=t\right)$ that depends on both the treatment assignment and the definition of the exposure mapping. Consequently, any change in the treatment assignment alters the estimand. As a result, ${\tau} = G \mu$ may lack a causal interpretation. One scenario in which the estimand $\tau$ can still be interpreted causally is when treatments are assigned independently and $T_i$ depends only on unit $i$'s local group of neighbors LeungLoupos2023. In this case, \( \tau \) represents a weighted average of unit-level treatment or spillover effects, comparing outcomes across different treatment assignments within this local group. For a more general discussion of causal inference with misspecified exposure mappings, see Savje2023.
To conclude this subsection, we present three examples of exposure mappings and interpret the corresponding estimands with some choices of $G$, in the context of PaluckShepherdAronow2016, which we will revisit in Section (ref). PaluckShepherdAronow2016 conducted a randomized experiment to study how an anti-conflict intervention influences teenagers' social norms regarding hostile behaviors such as bullying, social exclusion, harassment, and rumor-spreading. The treatment indicator $D_i$ corresponds to whether student $i$ was randomly assigned to participate in bi-weekly meetings that incorporated an anti-conflict curriculum. The outcome $Y_i$ is self-reported data on wristband wearing--a public signal of anti-conflict behavior and participation in the program. The network $A$ is measured by asking students to name up to ten students at the school they spent time with in the last few weeks.
Different specifications of the exposure mapping may change the estimand. For instance, the estimand defined using $T_i = T_{1i}$ or $T_i = T_{2i}$ alone differs from that obtained using a two-dimensional exposure mapping, $T_i = (T_{1i}, T_{2i})$, unless $T_{1i}$ and $T_{2i}$ are orthogonal. With independent $D_i$'s, the components $T_{1i}$ and $T_{2i}$ of the exposure mapping in Example (ref) are orthogonal. Therefore, the exposure mappings in Examples (ref) and (ref) respectively capture the direct and spillover effects in Example (ref). We examine all exposure mappings from Examples (ref)–(ref) when revisiting the empirical applications of PaluckShepherdAronow2016 and CaiJanvrySadoulet2015 to assess the robustness of our results to variations in the number of exposures; see Section (ref) and Appendix (ref).
Inverse probability weighting is a general estimation strategy in survey sampling and causal inference. In the context of observational studies with interference, TchetgenTchetgenVanderWeele2012, LiuHudgensBecker-Dreps2016 and JacksonLinYu2020 studied inverse probability-weighted estimators of causal effects under different assumptions on the interference pattern. In this subsection, we will review the Horvitz--Thompson and Hájek estimators for estimating population parameters based on the observed data in network experiments.
The Horvitz--Thompson estimator is a weighted estimator that assigns each unit a weight equal to the inverse of its selection probability. Recall $T_i = T(i, {D}, {A})$, and define the generalized propensity score Imbens2000 as $\pi_i(t) = \mathbb{P}( T_i=t )$. The value of the propensity score is known by design and can be determined through exact calculation or approximation using Monte Carlo AronowSamii2017. The Horvitz--Thompson estimator for $\mu(t)$ equals $\hat{Y}_{\text{ht}}(t) = n^{-1} \sum_{i=1}^n {1}(T_i=t) Y_i/\pi_i(t)$. The Horvitz--Thompson estimator is unbiased if the propensity score $\pi_i(t)$'s are non-zero and is consistent under additional regularity conditions. Leung2022 focused on $\tau(t,t')$ and examined the asymptotic properties of the Horvitz--Thompson estimator $\hat{\tau}_{\text{ht}}(t, t^{\prime}) = \hat{Y}_{\text{ht}}(t) - \hat{Y}_{\text{ht}}(t')$.
The Hájek estimator refines the Horvitz--Thompson estimator by normalizing the Horvitz--Thompson estimator by dividing it by the sum of the individual weights involved in its definition: $\hat{Y}_{\textup{haj}}(t) = \hat{Y}_{\text{ht}}(t)/\hat{1}_{\text{ht}}(t)$, where $\hat{1}_{\text{ht}}(t) = n^{-1} \sum_{i=1}^n {1}(T_i=t)/\pi_i(t)$ is the Horvitz--Thompson estimator for constant potential outcome $1$. The Hájek estimator is biased in general since $\hat{1}_{\text{ht}}(t)$ is random, but it is consistent since $\hat{1}_{\text{ht}}(t)$ is consistent for $1$ under regularity conditions.
The existing literature provides two motivations for using the Hájek estimator. First, it ensures invariance under the location shift of the outcome Fuller2011. Second, empirical evidence suggests that the Hájek estimator is more stable and efficient with little cost of bias in most reasonable scenarios LuncefordDavidian2004, Fuller2011, Ding2024. Leung2022 mentioned the Hájek estimator in the footnote of his paper without detailed theory. Moreover, the Hájek estimator is more natural from the regression perspective. Numerically, the Hájek estimator is identical to the coefficient from the WLS fit based on the inverse probability of the exposure mapping AronowSamii2017. The regression-based approach offers three notable advantages. First, WLS is easy to implement without too much additional programming. Second, WLS can provide network-robust standard errors. Third, WLS can incorporate covariates to improve efficiency when covariates are predictive of the outcome. The main focus of our paper is to explore the design-based properties of the Hájek estimators obtained through the regression-based method and associated HAC covariance estimator.
We consider Leung2022's framework of ANI. ANI refers to a situation where treatments assigned to individuals who are farther away from the focal unit have a diminishing effect on the focal unit's response, although the effect is not necessarily zero.
In this subsection, we provide an overview of the key assumptions outlined in Leung2022, which serve as the foundation for our analysis. These conditions ensure the theoretical properties of the regression-based point and covariance estimators. For readers more interested in practical applications, they have the option to skip this subsection during their initial reading and focus on the procedures and properties presented in Sections (ref) and (ref).
Let $\ell_{ {A}}(i, j)$ denote the path distance between units $i$ and $j$ within network $A$, representing the length of the shortest path connecting them. The path distance refers to the smallest number of edges that must be crossed to journey from unit $i$ to unit $j$ within the network. Furthermore, $\ell_A(i, j)$ is defined as $\infty$ if $i\ne j$ and no path exists between units $i$ and $j$ and defined as $0$ if $i = j$. For a specific unit $i$, its $K$-neighborhood, denoted by $\mathcal{N}(i, K; {A})= \{j \in \mathcal{N}_n: \ell_{ {A}}(i, j) \leq K \}$, includes the set of units within network $A$ that are at most at a path distance of $K$ from unit $i$. Define ${d}_{{\mathcal{N}(i, K; {A})}}=(d_j: j \in \mathcal{N}(i, K; {A}))$ and $ {A}_{\mathcal{N}(i, K; {A})}=(A_{kl}: k, l \in \mathcal{N}(i, K; {A}))$ as the subvector of $ {d}$ and subnetwork of $ {A}$ on $\mathcal{N}(i, K; {A})$, respectively.
Assumption (ref) requires the interference pattern of interest to be local, implying that the exposure mapping indicators are weakly dependent. Specifically, ${1}(T_i=t) \perp \!\!\! \perp {1}(T_j=t)$ if $\ell_{ {A}}(i, j)>2 K$ for some $K$. For instance, $K=0$ for the exposure mapping in Example (ref) and $K=1$ for both in Examples (ref) and (ref). Assumption (ref) requires the generalized propensity scores to be uniformly bounded between 0 and 1. Assumption (ref) imposes uniform boundedness on the potential outcomes.
Let ${D}^{\prime}$ be an IID copy of $ {D}$. Define ${D}^{(i, s)}= ({D}_{\mathcal{N}(i, s; {A})}, {D}_{\mathcal{N}_n \backslash \mathcal{N}(i, s; {A})}^{\prime})$ as the concatenation of the subvector of $ {D}$ on $\mathcal{N}(i, s; {A})$ and the subvector of ${D}^{\prime}$ on $\mathcal{N}_n \backslash \mathcal{N}(i, s; {A})$. Define
where the expectation is over the randomness of $D$ and $D'$ with all potential outcomes fixed. The interference, caused by distant individuals with a distance of more than $s$ from the subject, is measured as the largest expected change in any individual's potential outcome when altering the treatment assignments of those distant individuals. Mathematically, ANI assumes that as the distance $s$ approaches infinity, the largest value of $\theta_{n, s}$, taken over all feasible networks, converges to zero, which is formalized in Assumption (ref) below.
In simpler terms, Assumption (ref) stipulates that interference from distant individuals should vanish as the distance becomes large. We skip Assumption 5 in Leung2022, which is for showing consistency of the Horvitz--Thompson estimator, and proceed to Assumption (ref) below for the asymptotic normality of the Hájek estimator. Define
as the $k$-th moment of the $m$-neighborhood size within network $A$. For any $H, H' \subseteq \mathcal{N}_n$, define $\ell_A(H, H') = \min\{\ell_A(i, j) : i \in H, j \in H'\}$. Define
as the set of paired couples $(i, k)$ and $(j, l)$ such that the units within each couple are at most path distance $m$ apart from each other, and the two pairs are exactly path distance $s$ apart. Similarly, define
as the set of paired couples $(i, k)$ and $(j, l)$ such that the units within each couple are at most path distance $m$ apart from each other, and $i$ and $j$ are exactly path distance $s$ apart. In Assumption (ref), we replace $\sigma_n^2$ from Leung2022 with the matrix ${\Sigma}_\textup{haj}$:
Theorem (ref) below will show that ${\Sigma}_\textup{haj}$ defined in (ref) is the asymptotic covariance of the Hájek estimator of $\mu$. Based on the definition of $\theta_{n, s}$ in (ref) and Leung2022, we define
where $K$ is the constant from Assumption (ref) and $\lfloor s \rfloor$ is $s$ rounded down to the nearest integer. Assumptions (ref) and (ref) both posit that interference diminishes with path distance. Additionally, Assumption (ref) imposes further that for some sequence $m_n$, $\tilde{\theta}_{n, s}$ diminishes to zero at a sufficiently rapid rate relative to the size of the $m_n$-neighborhood, moreover, constraints on the growth of $m_n$-neighborhoods, and ensures that $\tilde{\theta}_{n, m_n}$ decays at an adequately fast pace. Moreover, Assumption (ref) is closely related to the conditions proposed in ChandrasekharJacksonMcCormick2024 to achieve the asymptotic normality of sums of dependent random variables. The three components of Assumption (ref) below are analogous to their Assumptions 1--3.
Assumption (ref) corresponds to Assumption 3.4 of KojevnikovMarmerSong2019, which limits the extent of dependence across units of $1(T_i=t)\pi_i(t)^{-1}(Y_i-\mu(t))$'s through restrictions on the network. Leung2022a verifies Assumption (ref) for networks with polynomial or exponential neighborhood growth rates. We impose Assumption (ref) to ensure the asymptotic normality of the Hájek estimator of $\mu$. We defer Assumption (ref), which ensures the consistency of covariance estimation, to Section (ref).
Let $z_i = ({1}(T_i=t): t\in\mathcal{T})$ be the vector of exposure mapping indicators. Motivated by the inverse probability weighting in the Hájek estimator, we consider the WLS fit:
Let $\hat{\beta}_{\textup{haj}}$ denote the estimtors of coefficients for $z_i$ in (ref). Define the concatenated Hájek estimator vector as $\hat{Y}_{\textup{haj}} = (\hat{Y}_{\textup{haj}}(t): t\in\mathcal{T})$. The numerical equivalence $\hat{\beta}_{\textup{haj}} = \hat{Y}_{\textup{haj}}$ is a well known result and shows the utility of WLS in reproducing the Hájek estimators AronowSamii2017, Ding2024. Theorem (ref) below states the asymptotic normality of $\hat{\beta}_{\textup{haj}}$.
Theorem (ref) ensures the consistency of $\hat{\beta}_{\textup{haj}}$ for estimating $\mu$ and establishes ${\Sigma}_{\textup{haj}}$ as the asymptotic sampling covariance of $\sqrt{n} (\hat{\beta}_{\textup{haj}} - {\mu})$.
The regression-based approach provides an estimator for the standard error via the same WLS fit. Denote the design matrix of the WLS fit in (ref) by an $n\times |\mathcal{T}|$ matrix ${Z} = (z_1,\ldots,z_n)^\top$, where its rows are the vectors $z_i$ for each unit $i\in\mathcal{N}_n$. Construct the weight matrix ${W}=\text{diag}\{ {w}_i: i = 1,\ldots, n\}$ by placing the weights $w_i$ along the diagonal. Let $ {Y} = (Y_1,\ldots, Y_n)$ denote the vector of the observed outcomes. Diagonalize the residual $e_i$'s from the same WLS fit to form the matrix $ {e}_{\textup{haj}} = \text{diag}\{e_i: i=1,\ldots,n\}$. Define
as the network-robust covariance estimator of $\hat{\beta}_{\textup{haj}}$, where $ {K}_n$ is a uniform kernel matrix with $(i,j)$th entry $K_{n,ij}={1}(\ell_{ {A}}(i, j)\le b_n)$. Here, choosing $b_n > 0$ places nonzero weight on pairs at most path distance $b_n$ apart from each other in the network $A$, which accounts for the network correlation. While (ref) adopts the form of an HAC estimator commonly used in spatial econometrics literature, our paper first discusses its design-based properties under the regression-based analysis for network experiments.
We follow the discussion in Leung2022 regarding the choice of the bandwidth $b_n$. Define the average path length, $\mathcal{L}({A})$, as the average value of $\ell_{{A}}(i, j)$ over all pairs in the largest component of ${A}$. Here, a component of a network refers to a connected subnetwork where all units within the subnetwork are disconnected from those outside of it. Let $\delta( {A})=$ $n^{-1} \sum_{i=1}^n \sum_{j=1}^n A_{i j}$ be the average degree. Leung2022 suggests choosing the bandwidth $b_n$ as follows:
where $\lfloor\cdot\rceil$ means rounding to the nearest integer. The choice of bandwidth $b_n$ is based on the following two reasons. First, $b_n$ is set to be at least equal to $2K$ to account for the correlation in $\{1(T_i=t)\}_{i=1}^n$ as per Assumption (ref). If the exposure mapping is correctly specified, we can simply choose $b_n=2K$. Second, (ref) chooses a bandwidth of logarithmic or polynomial order depending on the growth rates of the average $K$-neighborhood size. The logarithmic order in $b_n$ applies when the growth rate is approximately exponential in $K$ and polynomial order applies when the growth rate is approximately polynomial in $K$. Furthermore, Leung2022 justifies that the bandwidth in (ref) satisfies Assumption (ref)(b)–(d) under polynomial and exponential neighborhood growth rates. Since $K$ is researcher-defined, and $\mathcal{L}({A})$ and $\delta( {A})$ can be computed from the observed network data, $b_n$ in (ref) can be determined accordingly. To align with Leung2022, we also recommend that researchers report results for multiple bandwidths in a neighborhood of (ref) as a robustness check. We use the empirical application in Section (ref) as an illustrative example to demonstrate how to select the bandwidth.
We impose Assumption (ref), as introduced in Leung2022, to ensure the consistency of the covariance estimator, where $b_n$ is the bandwidth defined in (ref). Denote by \[ \mathcal{N}^{\partial}(i, s; {A})=\{j \in \mathcal{N}_n: \ell_{ {A}}(i, j)=s\} \] the $s$-neighborhood boundary of unit $i$, which is the set of units exactly at a distance of $s$ from $i$, and \[ M_n^{\partial}(s)=n^{-1} \sum_{i=1}^n |\mathcal{N}^{\partial}(i, s; {A})|, \] its average size across units.
Assumption (ref)(a) demonstrates the trade-off between restrictions on the network topology through $M_n^{\partial}(s)$ and the degree of interference through $\tilde{\theta}_{n, s}$. Assumption (ref)(b) and (d) regulate the bandwidth $b_n$ by imposing conditions on the first and second moments of the $b_n$-neighborhood size within network $A$. Assumption (ref)(d) is used to derive the asymptotic bias, which closely mirrors Assumption (ref) with $b_n$ and $\mathcal{J}_n(s, \cdot)$ in place of $m_n$ and $\mathcal{H}_n(s, \cdot)$, respectively. Assumption (ref) strongly depends on the structure of the underlying network. Leung2022 uses a mixture of formal and heuristic arguments to show that the bandwidth $b_n$ in (ref) satisfies Assumption (ref)(b)--(d) for networks with polynomial or exponential neighborhood growth rates.
Define ${\Delta}_{\textup{haj}}$ as an $n \times |\mathcal{T}|$ matrix with $(i,t)$th element ${\Delta}_{\textup{haj},it} = 1(T_i=t)\pi_i(t)^{-1}(Y_i - \mu(t)) - (\mu_i(t) - \mu(t) )$, and $M$ as an $n \times |\mathcal{T}|$ matrix with $(i,t)$th element $M_{it} = \mu_i(t)-\mu(t)$. Of interest is how this regression-based covariance estimator approximates the true sampling covariance from the design-based perspective.
We use $_*$ to indicate that ${\Sigma}_{*, \textup{haj}}$ is the “oracle” version of covariance estimator, which takes the form of a HAC estimator. Theorem (ref) first demonstrates that ${\Sigma}_{*, \textup{haj}}$ closely approximates the asymptotic covariance ${\Sigma}_{\textup{haj}}$ and then presents the asymptotic bias of the network-robust covariance estimator in estimating ${\Sigma}_{*, \textup{haj}}$. The bias term $R_{\textup{haj}}$ adopts the form of an HAC covariance estimator of the individual-level expected response. The covariance estimation is asymptotically exact with constant individual-level expected response under any exposure mapping value $t\in\mathcal{T}$, which is similar to the canonical results of Neyman1923 without interference. In some cases, the uniform kernel used in the network-robust covariance estimator $\hat{{V}}_{\textup{haj}}$ may not be positive semi-definite. This issue can result in an anti-conservative covariance estimator, which can in turn affect the accuracy of hypothesis testing and confidence intervals. We will address this issue in the next subsection. Now we end this subsection with a remark on the literature of HAC covariance estimators for network and spatial data.
There are four main concerns regarding the properties of the HAC variance estimator. First, it should ideally be non-negative in finite-samples, despite the kernel not always being positive semi-definite. Second, the HAC estimator is biased in a design-based setting, and it is desirable for the bias term to be asymptotically non-negative to ensure conservative inference. Third, HAC estimators often yield values that are too small in finite-samples compared with the true variance, leading to false discoveries. Finally, a computationally feasible bandwidth sequence is necessary for ensuring the consistency of the HAC estimator.
In this subsection, we tackle these issues by proposing a modification to the uniform kernel. Our proposed modification preserves the network-robustness of the covariance estimator while ensuring that it remains positive semi-definite and conservative. Let $Q_n \Lambda_n Q_n^{\top}$ be the eigendecomposition of $ {K}_n$. As ${K}_n$ is symmetric, all its eigenvalues are real. We define the adjusted kernel matrix by truncating the negative eigenvalues at $0$ as ${K}_n^{+}:= Q_n \max\{ \Lambda_n, 0\} Q_n^{\top}$, where the maximum is taken element-wise. Letting ${K}_n^{-} := Q_n |\min \{ \Lambda_n, 0\}| Q_n^{\top}$ with the minimum taken element-wise, we can also write $ {K}_n^{+} = {K}_n + {K}_n^{-}$. By construction, the matrix $ {K}_n^{\diamond}$ ($\diamond=+,-$) is positive semi-definite, and we denote the $(i,j)$th entry of $ {K}_n^{\diamond}$ as ${K}_{n,ij}^{\diamond}$. If \( K_n \) were positive semi-definite, then $K_n = {K}_n^+ $. We propose the adjusted HAC covariance estimator as
To guarantee the asymptotic conservativeness of $\hat{ {V}}_{\textup{haj}}^+$, we impose Assumption (ref) below, which pertains to the properties of $K_n^{-}$. Recall that $K_{n,ij} = {1}(\ell_{ {A}}(i, j)\le b_n)$ and write $M_n(m,k)$ and $\mathcal{J}_n(s, m)$ in (ref) and (ref) with $m=b_n$ as:
Define $M_n^{-}(b_n,k)$ and $\mathcal{J}^{-}_n(s, b_n)$ as the counterparts of $M_n(b_n,k)$ and $\mathcal{J}_n(s, b_n)$ on $|K_n^-|$, respectively:
Assumption (ref) is the analogue of Assumption (ref), but specifically tailored to the quantity $|K_n^-|$, with Assumption (ref)(a) identical to Assumption (ref)(a).
Theorem (ref) delineates two key advantages stemming from the construction of the adjusted covariance estimator. First, it ensures that the covariance estimator $\hat{ {V}}_{\textup{haj}}^+$ is positive definite. Second, it produces a positively adjusted bias term $R_{\textup{haj}}^+$, leading to the conservativeness of $\hat{ {V}}_{\textup{haj}}^+$ for estimating the true sampling covariance. Theorems (ref) and (ref) together justify the regression-based inference of $\tau = G\mu$ from the WLS fit (ref) with the point estimator $\hat{\tau} = G\hat{\beta}_{\textup{haj}}$ and the adjusted regression-based HAC covariance estimator $G \hat{ {V}}_{\textup{haj}}^+ G^\top$.
In this subsection, we briefly discuss other covariance estimation strategies. KojevnikovMarmerSong2021 provides a law of large numbers and a central limit theorem for network dependent variables. Additionally, they introduce a technique for computing standard errors that remains robust under various types of network dependencies. Their approach relies on a network HAC covariance estimator for a broad class of kernel functions, which they show consistently estimates the true sampling covariance. As demonstrated in Leung2022, the uniform kernel provides better size control, especially in cases with smaller samples, compared with alternative kernels that diminish with distance. Considering these reasons, we opt for the uniform kernel.
Leung2022 proposes a covariance estimator for the Horvitz--Thompson estimator of exposure effects, while Kojevnikov2021 develops bootstrap-based alternatives to network HAC estimation. Although both estimators share similarities with the HAC framework, neither is derived from a regression-based approach. Kojevnikov2021 ensures that the resulting estimator is positive semi-definite, and Leung2019e refines this approach by showing that, under a specific bandwidth choice, the variance estimator exhibits non-negative asymptotic bias. However, both methods suffer from substantial overrejection in finite-sample simulations. We compare the finite-sample performance of our estimator with those of Leung2019e and Kojevnikov2021 in Section (ref).
Leung2022 and our regression-based HAC estimator \( \hat{V}_{\textup{haj}} \) both use the uniform kernel, which helps mitigate overrejection in finite-samples. As shown in Leung2022a, Leung2022's variance estimator is asymptotically conservative under mild weak dependence conditions on the super-population. However, the non-positive semi-definiteness of the uniform kernel can lead both estimators to produce negative variance estimates in finite-samples, resulting in potential anti-conservativeness in both asymptotic theory and simulations. The idea of replacing the negative eigenvalues of $K_n$ with non-negative values appeared in Kojevnikov2021, which can be traced back to the literature on approximating a symmetric matrix by a positive definite matrix Higham1988,Politis2009. The key distinction is that Kojevnikov2021 applied this technique to the final HAC covariance estimator, while we apply it to the kernel matrix. There are two limitations of Kojevnikov2021's approach. First, it is not suitable for estimating a single causal effect, as when the HAC estimator is scalar, it merely involves replacing a negative variance estimate with zero. In contrast, our approach is applicable to joint causal effects. Second, Kojevnikov2021's approach does not address the issue of anti-conservativeness, as the crucial factor for positive bias is the positive semi-definiteness of $K_n$. WangSamiiChang2025 recently applied our strategy to the HAC variance estimator in the spatial experiments and found better finite-sample properties.
Regression-based methods offer a natural framework for incorporating covariates and can lead to efficiency gains under appropriate conditions.\footnote{AronowSamii2017 discussed the use of covariates to improve efficiency via difference estimators, although they did not implement this approach in their analysis.} To set the stage for our discussion, we briefly review the theory of covariate adjustment under complete randomization without interference.
Consider an experimental setup involving a binary intervention and a population of $n$ units with potential outcomes denoted by $Y_i(0)$ and $Y_i(1)$ for each unit $i = 1, \ldots, n$. The average treatment effect within the finite population is denoted by $\tau(1,0) = \bar{Y}(1) - \bar{Y}(0)$, where $\bar{Y}(z) = n^{-1}\sum_{i=1}^n Y_i(z)$ for $z = 0, 1$. Denote by $z_i$ the treatment indicator of unit $i$ under complete randomization. The difference-in-means estimator is unbiased for $\tau(1,0)$, and equals the coefficient of $z_i$ from the Ordinary Least Squares (OLS) regression of $Y_i$ on $(1, z_i)$. Given the covariate vector $x_i = (x_{i1}, \ldots, x_{iJ})$ for $i = 1, \ldots, n$, Fisher1935 proposed to use the coefficient of $z_i$ from the OLS fit of regressing $Y_i$ on $(1,z_i,x_i)$ to estimate $\tau(1,0)$. Freedman2008 criticized this approach, highlighting its potential for efficiency loss compared to the difference-in-means estimator. Lin2013 introduced an improved estimator, defined as the coefficient of $z_i$ obtained from the OLS regressing of $Y_i$ on $(1,z_i,(x_i-\bar{x}), z_i(x_i-\bar{x}))$. This specification includes covariates as well as treatment-covariate interactions. He proved that this estimator is at least as efficient as the difference-in-means and Fisher1935's estimators in the asymptotic sense.
We refer to the regression proposed by Fisher1935 as the additive specification, and Lin2013's regression as the fully-interacted specification to avoid any ambiguity. We expand upon their findings in the context of network experiments, which incorporate interference, through the utilization of WLS fits. To simplify the presentation, we center the covariates at $\bar{x} = n^{-1} \sum_{i=1}^n x_i = 0$.
Recall $z_i = (1(T_i=t): t\in \mathcal{T})$ as the dummies for the exposure mapping in the network experiment. Consider the WLS fit
Let $\hat{\beta}_{\textup{haj},\textsc{f}}$ denote the estimtors of coefficients for $z_i$ from the above WLS fit and $\hat{\beta}_{\textup{haj},\textsc{f}}(t)$ denote the element in $\hat{\beta}_{\textup{haj},\textsc{f}}$ corresponding to $1(T_i=t)$. We use the subscript “F” to signify Fisher1935. Assumption (ref) below imposes the uniform boundedness of $x_i$ and adapts Assumption (ref) to its version with covariate adjustment.
Let ${\gamma}_\textsc{f}$ denote the probability limit of $\hat{\gamma}_\textsc{f}$, where $\hat{\gamma}_\textsc{f}$ is the coefficient vector of $x_i$ from the WLS fit in (ref). Let ${\Sigma}_{\textup{haj},\textsc{f}}$ denote the analog of ${\Sigma}_{\textup{haj}}$ in (ref) defined on the covariate-adjusted outcome $Y_i-x_i^\top {\gamma}_\textsc{f} $. Theorem (ref) below states the asymptotic normality of $\hat{\beta}_{\textup{haj},\textsc{f}}$.
The design matrix of the WLS fit in (ref) equals $ {C}_\textsc{f} = ( {Z}, {X})$ where $Z$ is an $n\times |\mathcal{T}|$ matrix and $ {X} = (x_i:i=1,\ldots,n)$ is an $n \times J$ matrix. Diagonalize the residual ${e_{\textsc{f},i}}$'s from the WLS fit in (ref) to form the matrix $ {e}_{\textup{haj},\textsc{f}} = \text{diag}\{e_{\textsc{f},i}: i=1,\ldots, n\}$. Let $[\cdot]_{(1:|\mathcal{T}|,1:|\mathcal{T}|)}$ denote the upper-left $|\mathcal{T}|\times |\mathcal{T}|$ submatrix. Let $\hat{{V}}_{\textup{haj},\textsc{f}}$ denote the HAC estimator for $\hat{\beta}_{\textup{haj},\textsc{f}}$, which is a submatrix of the covariance estimator obtained from the WLS fit in (ref): \[ \hat{ {V}}_{\text{haj,\textsc{f}}} = \left[( {C}_\textsc{f}^{\top} {W} {C}_\textsc{f})^{-1} ( {C}_\textsc{f}^{\top} {W} {e}_{\textup{haj},\textsc{f}} {K}_n {e}_{\textup{haj},\textsc{f}} {W} {C}_\textsc{f}) ( {C}_\textsc{f}^{\top} {W} {C}_\textsc{f})^{-1} \right]_{(1:|\mathcal{T}|,1:|\mathcal{T}|)}. \] Let ${\Delta}_{\textup{haj}, \textsc{f}}$ denote the analog of ${\Delta}_{\textup{haj}}$ defined on the covariate-adjusted outcome $Y_i-x_i^\top {\gamma}_\textsc{f} $. Define $M_\textsc{f}$ as an $n\times |\mathcal{T}|$ matrix with $(i,t)$th element $M_{\textsc{f}, it} = \mu_i(t)-\mu(t) - x_i^\top \gamma_\textsc{f}$. Theorem (ref) below establishes the asymptotic bias of $\hat{{V}}_{\text{haj,\textsc{f}}}$ as an estimator for the asymptotic covariance of $\hat{\beta}_{\text{haj,\textsc{f}}}$.
The bias term $R_{\textup{haj},\textsc{f}}$ is an analog of $R_{\textup{haj}}$ defined on the adjusted outcome $Y_i-x_i^\top {\gamma}_\textsc{f} $. Given that $K_n$ may not be positive semi-definite, we cannot ensure the asymptotic conservativeness of $\hat{ {V}}_{{\textup{haj},\textsc{f}}}$ for estimating ${ {\Sigma}}_{*,\textup{haj}, \textsc{f}}$. Similar to (ref), we propose the adjusted covariance estimator as \[ \hat{ {V}}_{\textup{haj},\textsc{f}}^+ = \left[( {C}_\textsc{f}^{\top} {W} {C}_\textsc{f})^{-1} ( {C}_\textsc{f}^{\top} {W} {e}_{\textup{haj},\textsc{f}} {K}_n^+ {e}_{\textup{haj},\textsc{f}} {W} {C}_\textsc{f}) ( {C}_\textsc{f}^{\top} {W} {C}_\textsc{f})^{-1} \right]_{(1:|\mathcal{T}|,1:|\mathcal{T}|)}. \]
Theorem (ref) ensures the asymptotic conservativeness of $\hat{ {V}}_{\text{haj,\textsc{f}}}^+$ for estimating the true sampling covariance. This, together with Theorem (ref), justify the regression-based inference of $\tau = G\mu$ from the additive WLS fit in (ref) with the point estimator $\hat{\tau} = G\hat{\beta}_{\text{haj,\textsc{f}}}$ and the adjusted regression-based HAC covariance estimator $G \hat{ {V}}_{\text{haj,\textsc{f}}}^+ G^\top$.
With full interactions between the exposure mapping indicators and covariates, we consider the WLS fit
where $\otimes$ denotes the Kronecker product. The specification (ref) simply means WLS fit of $Y_i$ on the dummy $1(T_i=t)$'s and the interaction $1(T_i=t)x_i$'s. Let $\hat{\beta}_{\textup{haj},\textsc{l}}$ denote the estimtors of coefficients for $z_i$ from the above WLS fit and $\hat{\beta}_{\textup{haj},\textsc{l}}(t)$ denote the element in $\hat{\beta}_{\textup{haj},\textsc{l}}$ corresponding to $1(T_i=t)$. We use the subscript “L” to signify Lin2013. Let ${\gamma}_\textsc{l}(t)$ be the probability limit of $\hat{\gamma}_\textsc{l}(t)$, where $\hat{\gamma}_\textsc{l}(t)$ is the coefficient vector of $1(T_i=t)x_i$ from the WLS fit in (ref). Let ${\Sigma}_{\textup{haj},\textsc{l}}$ be the analog of ${\Sigma}_{\textup{haj}}$ in (ref) defined on the adjusted outcome $Y_i- x_i^\top {\gamma}_\textsc{l}(T_i)$. Theorem (ref) states the asymptotic normality of $\hat{\beta}_{\textup{haj},\textsc{l}}$.
Let $ {C}_\textsc{l}$ be the design matrix of the WLS fit in (ref), with row vectors $(z_i^\top, (z_i \otimes x_i)^\top)$. Diagonalize the residual $e_{\textsc{l},i}$'s from the same WLS fit to form the matrix $ {e}_{\textup{haj},\textsc{l}} = \text{diag}\{e_{\textsc{l},i}:i=1,\ldots,n\}$. Let $\hat{ {V}}_{\textup{haj},\textsc{l}}$ denote the HAC covariance estimator for $\hat{\beta}_{\textup{haj},\textsc{l}}$, which is a submatrix of the covariance estimator obtained from the WLS fit in (ref): \[ \hat{ {V}}_{\text{haj,\textsc{l}}} = \left[( {C}_\textsc{l}^{\top} {W} {C}_\textsc{l})^{-1} ( {C}_\textsc{l}^{\top} {W} {e}_{\textup{haj},\textsc{l}} {K}_n {e}_{\textup{haj},\textsc{l}} {W} {C}_\textsc{l}) ( {C}_\textsc{l}^{\top} {W} {C}_\textsc{l})^{-1} \right]_{(1:|\mathcal{T}|,1:|\mathcal{T}|)}. \] Let $ {\Delta}_{\textup{haj}, \textsc{l}}$ be the analog of $ {\Delta}_{\textup{haj}}$ defined on the adjusted outcome $Y_i- x_i^\top {\gamma}_\textsc{l}(T_i)$. Define $M_\textsc{l}$ as an $n\times |\mathcal{T}|$ matrix with $(i,t)$th element $M_{\textsc{l}, it} = \mu_i(t)-\mu(t) - x_i^\top \gamma_\textsc{l}(t)$.
Theorem (ref) establishes the asymptotic bias of $\hat{ {V}}_{{\textup{haj},\textsc{l}}}$ as an estimator for the asymptotic covariance of $\hat{\beta}_{\text{haj,\textsc{f}}}$. Given that $K_n$ may not be positive semi-definite, we cannot ensure the asymptotic conservativeness of $\hat{ {V}}_{{\textup{haj},\textsc{l}}}$ for estimating $\hat{ {\Sigma}}_{*,\textup{haj}, \textsc{l}}$. Similar to (ref), we propose the adjusted HAC covariance estimator as \[ \hat{ {V}}_{\textup{haj},\textsc{l}}^+ = \left[( {C}_\textsc{l}^{\top} {W} {C}_\textsc{l})^{-1} ( {C}_\textsc{l}^{\top} {W} {e}_{\textup{haj},\textsc{l}} {K}_n^+ {e}_{\textup{haj},\textsc{l}} {W} {C}_\textsc{l}) ( {C}_\textsc{l}^{\top} {W} {C}_\textsc{l})^{-1} \right]_{(1:|\mathcal{T}|,1:|\mathcal{T}|)}. \]
Echoing the comment after Theorem (ref), Theorems (ref) and (ref) together justify the regression-based inference of $\tau = G\mu$ from the fully-interacted WLS fit in (ref) with point estimator $\hat{\tau} = G\hat{\beta}_{\text{haj,\textsc{l}}}$ and adjusted regression-based HAC covariance estimator $G \hat{ {V}}_{\text{haj,\textsc{l}}}^+ G^\top$.
Regression adjustment can improve efficiency under reasonable data-generating processes. Lin2013 demonstrated the efficiency gain from including fully interacted covariates when the propensity score is constant and there is no interference. However, this strategy does not always improve efficiency, especially in the presence of heterogeneous propensity scores or interference. In Appendix (ref), we present simulation results demonstrating that including fully-interacted covariates can exhibit higher asymptotic variance than the unadjusted Hájek estimator in scenarios with either heterogeneous propensity scores or interference. This lack of guarantee has also been documented in settings without interference, such as cluster experiments with varying sizes SuDing2021, split-plot experiments ZhaoDing2022, and scenarios where outcomes are not missing completely at random ZhaoDingLi2024. Despite the lack of theoretical guarantees for efficiency gain, we do observe that covariate adjustment improves efficiency in the simulation studies and empirical examples in Section (ref).
We focus on regression-based covariate-adjusted estimators for ease of implementation. There are alternative methods for enhancing efficiency via covariate adjustment. One strategy is to find the optimal linearly adjusted estimator by minimizing the true or estimated standard error; see, e.g., LiDing2020a and LuShiFang2025. For example, in our setting, we consider regressions such as regressing \( Y_i - x_i^\top \gamma \text{ on } z_i \) or regressing \( Y_i - 1(T_i = t) x_i^\top \gamma(t) \text{ on } 1(T_i = t) \) for each $t\in \mathcal{T}$. We can compute the variance of the resulting estimator and then minimize it with respect to the coefficients $\gamma$ or \( \gamma(t) \)'s. This procedure is different from WLS, but the minimization ensures variance reduction. Another strategy is to first compute the Hájek estimators of \( Y_i \) and \( x_i \), denoted by \( \hat{\beta}_{\textup{haj}} \) and \( \hat{\beta}_{{\textup{haj}, x}} \), respectively. An adjusted estimator can then be constructed as $\hat\beta_{\textup{adj}} = \hat\beta_{\textup{haj}} - \hat\beta_{{\textup{haj}, x}}^\top \gamma_x$, where the optimal adjustment coefficient is given by $\gamma_x = \operatorname{Cov}(\hat{\beta}_{{\textup{haj}, x}})^{-1} \operatorname{Cov}(\hat{\beta}_{{\textup{haj}, x}}, \hat\beta_{\textup{haj}})$, and can be consistently estimated; see JiangandWei2019 and RothSantAnna2023. When using the estimated covariance matrix of $\hat{\beta}_{\textup{haj}}$ and $\hat\beta_{{\textup{haj}, x}}$ to estimate $\gamma_{x}$, the resulting estimated standard error is always smaller. We omit the details for these two alternative strategies because we focus on simpler regression-based estimators.
In this section, we first examine the finite-sample performance of our results with simulation and then apply our results to an empirical application. Our analysis focuses on the exposure effect $\tau(t,t') = \mu(t) - \mu(t')$. We analyze another empirical example CaiJanvrySadoulet2015 in Appendix (ref).
To achieve comparability with Leung2022, we replicate the same scenario but with the inclusion of a covariate in the model. Regarding the results, we present the point and covariance estimators of the exposure effect from three specifications of WLS: unadjusted (Unadj), with additive covariates (Add), and with fully-interacted covariates (Sat). We also report Leung2022's Horvitz--Thompson estimator and variance estimator.
The study encompasses two outcome models: the linear-in-means model and the complex contagion model. Define
where $\tilde{A}_{i j} = A_{ij}/\sum_{j=1}^n A_{i j}$ is the $(i,j)$th entry of $\tilde{A}$, the row-normalized version of $A$. For the linear-in-means model, we set $Y_i = V_i( {D}, {A}, {x}, {\varepsilon})$ with $(\alpha, \beta, \delta, \xi, \gamma)=(-1,0.8,1,1,3)$. The model defines potential outcomes $Y_i(D)$ through its reduced form: \[ Y = \alpha (I - \beta \tilde{A})^{-1} \iota + (I - \beta \tilde{A})^{-1} (\delta \tilde{A} + \xi I) D + (I - \beta \tilde{A})^{-1} \gamma x + (I - \beta \tilde{A})^{-1} \varepsilon. \] For the complex contagion model, we set $Y_i={1}(V_i( {D}, {A}, {x}, {\varepsilon})>0)$ with $(\alpha, \beta, \delta,\xi, \gamma)=(-1,1.5,1,1,3)$. The complex contagion model can be generated from the dynamic process: \[ Y_i^t = 1\left( \alpha+\beta \sum_{j=1}^n \tilde{A}_{i j} Y_j^{t-1} + \delta \sum_{j=1}^n \tilde{A}_{i j} D_j + \xi D_i + \gamma x_i + \varepsilon_i > 0 \right) \] with initialization at period 0 as
We run the dynamic process to obtain new outcomes $Y^t=(Y_i^t)_{i=1}^n$ from last period's outcomes $Y^{t-1}$ until the first period $T$ such that $Y^T = Y^{T-1}$. We then take $Y^T$ as the vector of observed outcomes $Y$, which yields outcomes $(Y_i(D))^n_{i=1}$. As a result, this process implicitly defines potential outcomes (Leung2022). Without covariates $x_i$, Leung2022 derived conditions on the model parameters of the linear-in-means model and complex contagion model so that ANI holds. We can extend his proof to the models with additive covariates as in (ref), or covariates interacted with the network $A$, given that the covariates are fixed. We choose parameters to satisfy those conditions to ensure ANI.
Following Leung2022, we generate the adjacency matrix $ {A}$ from a random geometric graph model. Specifically, for each node $i$, we randomly generate its position $\rho_i$ in a two-dimensional space from $\mathcal{U}([0,1]^2)$. An edge between nodes $i$ and $j$ is created if the Euclidean distance between their positions is less than or equal to a threshold value $r_n$: $A_{i j}={1}\{\left\|\rho_i-\rho_j\right\| \leq r_n\}$, where the threshold value is chosen as $r_n=(\kappa /(\pi n))^{1/2}$. We set $\kappa$ as the average degree $\delta({A})$, calculated based on the experimental data in Section (ref), in order to better mimic real-world scenarios. We also generate a sequence $\{\nu_i\}_{i=1}^n \stackrel{\text{IID}}{\sim} \mathcal{N}(0,1)$ independent of $ {A}$. The error term in (ref) is generated as $\varepsilon_i=\nu_i + (\rho_{i 1}-0.5)$, where $\rho_{i 1}$ is the first component of $i$'s “location” $\rho_{i}$ generated above. This inclusion accounts for unobserved homophily, as units with similar $\rho_{i1}$ values are more likely to form links. Finally, we generate the covariate $\{x_i\}_{i=1}^n \stackrel{\text{IID}}{\sim} \mathcal{N}(0,1)$.
We use the sample of the two largest treated schools from the network experiment in Section (ref) to calibrate the network models. The network size $n$ is $1456$. We also conduct simulation with network sizes $n=805$ and $2725$ to illustrate variations in population sizes. See results in the Appendix (ref). We treat the schools as a single network by pooling the degree sequences across them. We randomly assign treatments to units classified as eligible in the experimental data with a probability $0.5$. Since we work within a finite-population framework, we generate $ {A}$, $ {\varepsilon}$'s, and $ {x}$'s once and only redraw $ {D}$ for each simulation draw. This differs from the superpopulation design simulation in Leung2022, where he regenerated $D$, $A$ and $\varepsilon$'s for each simulation draw.
For the spillover effect of having at least one treated friend versus non-treated friends $\tau(1,0)$, we define the exposure mapping as $T_i={1}(\sum_{j=1}^n A_{i j} D_j>0)$ and analyze only the population of units with at least one friend who is eligible for treatment to satisfy Assumption (ref). Under the IID randomization of $D$, we can compute the propensity score $\pi_i(1)$'s and $\pi_i(0)$'s for each student using Binomial probabilities.
Table (ref) presents the results. The top panels display our regression-based results. We report the estimand under “${\tau}(1,0)$,” approximated by the unbiased Horvitz--Thompson estimator $\hat{\tau}_\text{ht}(1,0)$, computed over $10,000$ simulation draws. We report “Oracle SE,” denoted by $\operatorname{Var}(\hat{\tau}(1,0))^{1/2}$, which are calculated as the standard deviation of the point estimators from corresponding WLS fits over $10,000$ simulation draws. For the estimation results, we conduct another independent $10,000$ simulation draws. We present the point estimate from each WLS fit under “$\hat{\tau}(1,0)$.” We present the HAC standard errors obtained from each WLS fit under “WLS SE,” and the corresponding adjusted HAC standard errors under “WLS$^+$ SE”, where the suggested bandwidth based on (ref) is $b_n = 3$. We report the Eicker-Huber-White standard errors assuming no interference under “EHW SE” to illustrate the degree of dependence in the data. We also report the empirical coverage rate of $95\%$ confidence intervals (CIs) in the “Coverage” rows for the corresponding standard errors. The effective sample size of exposure mapping value $t$ is defined as $\hat{n}(t)=\sum_{i=1}^n 1(T_i=t)$.
The result table demonstrates that the standard errors obtained from the WLS fits can be anti-conservative, underestimating the true standard error. However, by utilizing the adjusted HAC standard errors, we can improve the empirical coverage and ensure a conservative estimation of the standard error. In this setting, the estimator from the fully-interacted WLS fit is at least as efficient as the estimators from the unadjusted or additive WLS fits.
In the middle panel of Table (ref), we report the results of standard errors and coverage rates of $95\%$ CIs using the kernel ${K}^{\textup{L}2019}_n$ in Leung2019e and the kernel ${K}^{\textup{K}2021}_n$ in Kojevnikov2021. Both ${K}^{\textup{L}2019}_n$ and ${K}^{\textup{K}2021}_n$ are positive semi-definite, ensuring the positive semi-definiteness of the covariance estimators. However, we can see that they substantially overreject even in moderately sized samples.
The bottom panel of Table (ref) present the results of the Horvitz--Thompson estimator and variance estimator from Leung2022. By comparing the “Oracle SE” from the top and bottom panels, we can see the WLS estimators from all three specifications exhibit higher efficiency compared with the Horvitz--Thompson estimator. Moreover, Leung2022’s standard errors are smaller than the oracle standard errors, resulting in under coverage.
In this subsection, we revisit PaluckShepherdAronow2016 and apply our regression-based analysis to their network experiment, which examines how an anti-conflict intervention influences teenagers’ social norms regarding hostile behaviors such as bullying, social exclusion, harassment, and rumor-spreading. We now provide a detailed description of the empirical setting. In the experimental design, half of $56$ schools were randomly assigned to the treatment group. Within these treated schools, a subset of students was selected as eligible for treatment based on certain characteristics. Half of the eligible students were then block-randomized into treatment by gender and grade. Those treated students were invited to participate in bi-weekly meetings that incorporated an anti-conflict curriculum. Following Leung2022, we choose self-reported data on wristband wearing as the outcome of interest, which serves as the reward for students who exhibit anti-conflict behavior. We incorporate both gender and grade for covariate adjustment. The network is measured by asking students to name up to ten students at the school they spent time with in the last few weeks. More details about this network experiment can be found in PaluckShepherdAronow2016.
To align with the results reported in Leung2022, we restrict the data to the five largest treated schools. Our primary interest lies in assessing the direct effect of the anti-conflict intervention and the spillover effect of having at least one friend assigned to the treatment versus none such friends. We first calculate both effects by defining two one-dimensional exposure mappings and report the results in Table (ref). To examine both effects simultaneously, we define a two-dimensional exposure mapping and report the results in Table (ref). The network, obtained from surveys, is directed. When calculating the number of treated friends for the exposure mappings, we take into account the direction of links. However, when computing network neighborhoods for our covariance estimators, we disregard the directionality of links to conservatively define larger neighborhoods. For each exposure mapping, our analysis involves three WLS specifications: unadjusted (Unadj), with additive covariates (Add), and with fully-interacted covariates (Sat). We also include the results from Leung2022 in the column “Leung.”
\paragraph*{One-dimensional exposure mapping} For the direct effect, we define $T_i=D_i$ as in Example (ref) and limit the analysis to the students eligible for treatment, totaling $320$ students. The propensity score is $\pi_i(t)=0.5$ for each student. For the spillover effect, we employ $T_i={1}(\sum_{j=1}^n A_{i j} D_j>0)$ as the exposure mapping as in Example (ref), indicating whether at least one friend has been assigned to the treatment. We restrict the effective sample to units with at least one eligible friend. Under block randomization, we can compute the propensity score $\pi_i(0)$ and $\pi_i(1)$ for each student using Hypergeometric probabilities.
The results are presented in Table (ref). The suggested bandwidths in (ref) are $b_n = 2$ for both exposure mappings. We present results for the range of bandwidths $\{0, \ldots, 3\}$, where $0$ yields the standard errors in the absence of interference. The first row, labeled as “Estimate,” presents the point estimates obtained from corresponding WLS fits. The rows labeled as “$b_n=k$” present the HAC standard errors with the specific bandwidth values stated. We find that the kernel matrix $K_n$ is not positive semi-definite for all bandwidths in $\{1, 2, 3\}$, so we report the adjusted HAC standard errors under “WLS$^{+}$ SE”. The direct effect is statistically significant at $5\%$ level across all specifications, bandwidths, and after adjustment to the covariance estimation. The spillover effect is significant at $5\%$ level except when $b_n = 3$, both before and after adjustment to the covariance estimation. While our results align with the conclusions of Leung2022, our regression-based estimation approach provides higher precision. Also, the $K_n$ is not positive semi-definite indicating that Leung2022's variance estimators may be anti-conservative.
\paragraph*{Two-dimensional exposure mapping} We define the exposure mapping and $G$ as in Example (ref): $T_i = (D_i, {1} (\sum_{j=1}^n A_{ij}D_j > 0) )$. We focus on the first two components of $\tau=G\mu$, where the first component captures the direct effect and the second component captures the spillover effect. We restrict the effective sample to students who are eligible for treatment and have at least one eligible friend, resulting in a total of $150$ students.
The results are presented in the top panel of Table (ref). The average out-degree, $n^{-1} \sum_{ij} A_{ij}$, is 7.96. The APL is 3.37 across our five schools. Given $n = 3306$ students, we have $\log n / \log \delta(A) = 3.96$, which is close to 3.37. Thus, the suggested bandwidth in (ref) is $b_n = 2$ with $K = 1$, and we report results for the range of bandwidths $\{0, \ldots, 3\}$. We observe that the magnitude and standard errors of the direct effect remain relatively stable. Regarding the spillover effect, its magnitude notably increases, and it remains statistically significant at the $5\%$ significance level across all specifications and bandwidths, even after adjustment to the covariance estimation.
To investigate whether these changes in results arise from shifts in the target population or potential misspecification of the exposure mappings, we provide results using two one-dimensional exposure mappings and focusing on treatment-eligible students with at least one eligible friend. These results are displayed in the bottom panel of Table (ref). Upon comparing the top and bottom panels, we can observe that there are minor differences in the point estimates and standard errors, but the overall message does not change. Specifically, the spillover effect is more pronounced and significant for the subset of students who are both eligible for treatment and have at least one eligible friend, in comparison to the subset with at least one eligible friend. Table (ref) also demonstrates that our methods are robust to various specifications of exposure mappings.
Our theory focuses on discrete exposure mappings with finite support. However, continuous or growing-dimensional exposure mappings, such as the number or share of treated friends, are also common in practice, e.g., MuralidharanNiehausSukhtankar2023. For growing-dimensional exposure mappings that vary with $n$, valid inference is possible when the network is sparse, meaning that the maximum or average degree is substantially smaller than the network size (e.g., Leung2020). For continuous exposure mappings, estimating \( \mu(t) \) is conceptually straightforward by extending the propensity score to a treatment density function, defined as \( \pi_i(t) = f_{T_i}(t) \) HiranoImbens2004.
Without imposing any modeling assumption on $\mu(t)$, we can use the following nonparametric estimator:
which locally averages the $Y_i$ values whose $T_i$ falls within the bandwidth $h$ around $t$. Since the exposure mapping $T_i$ and the treatment assignments are known, one can compute $\mathbb{P}(|T_i - t| \leq h)$ either in closed form or via Monte Carlo simulation. The main technical challenges are (i) ensuring sufficient smoothness of the estimand $\mu(t)$ and (ii) choosing an appropriate bandwidth $h$ to trade off the bias and variance for estimating $\mu(t)$. Here, we use the uniform kernel in (ref) as an illustrative example, although general kernel functions could be employed.
As noted by FaridaniNiehaus2024, regression-based analysis with continuous exposure mappings typically relies on either a linear outcome model or restrictions on the experimental design. Consider the potential outcome model \( Y_i(t) = Y_i(0) + \beta_i t \), where the individual effects $\beta_i$'s can vary across units. If we regress the outcome \( Y_i \) on the centered exposure mapping \( T_i - \mathbb{E}(T_i) \) with weight $1/\sqrt{\operatorname{Var}(T_i)}$, the WLS coefficient
identifies: \[ \beta = \frac{\frac{1}{n} \sum_{i=1}^n \frac{1}{\operatorname{Var}(T_i)} \operatorname{Cov}(Y_i, T_i - \mathbb{E}(T_i))}{\frac{1}{n} \sum_{i=1}^n \frac{1}{\operatorname{Var}(T_i)} \operatorname{Var}(T_i)} = \frac{1}{n} \sum_{i=1}^n \beta_i, \] which represents the average of the $\beta_i$'s. With constant treatment effect $\beta_i = \beta$, the WLS coefficient $\hat{\beta}$ identifies $\beta$.
We outline future directions for continuous exposure mapping above, leaving many technical issues for further research. For example, what is the optimal choice of bandwidth $h$ in estimator (ref)? More importantly, we aim to develop rigorous statistical inference procedures for both the nonparametric estimator in (ref) and the WLS estimator in (ref).