# Almost-Matching-Exactly for Treatment Effect Estimation under Network
Interference

M. Usaid Awan  Department of Economics, Duke University  Marco Morucci
Department of Political Science, Duke University  Vittorio Orlandi  Department
of Statistical Science, Duke University  Sudeepa Roy  Department of Computer
Science, Duke University  Cynthia Rudin  Department of Statistical Science,
Duke University  Department of Computer Science, Duke University  Department
of Electrical and Computer Engineering, Duke University  Alexander Volfovsky
Department of Statistical Science, Duke University

###### Abstract

We propose a matching method that recovers direct treatment effects from
randomized experiments where units are connected in an observed network, and
units that share edges can potentially influence each others’ outcomes.
Traditional treatment effect estimators for randomized experiments are biased
and error prone in this setting. Our method matches units almost exactly on
counts of unique subgraphs within their neighborhood graphs. The matches that
we construct are interpretable and high-quality. Our method can be extended
easily to accommodate additional unit-level covariate information. We show
empirically that our method performs better than other existing methodologies
for this problem, while producing meaningful, interpretable results.

###  1 INTRODUCTION

Randomized experiments are considered to be the gold standard for estimating
causal effects of a treatment on an outcome. Typically, in these experiments,
the outcome of a unit is assumed to be only affected by the unit’s own
treatment status, and not by the treatment assignment of other units (Cox,
1958; Rubin, 1980). However, in many applications – such as measuring
effectiveness of an advertisement campaign or a teacher training program –
units interact, and ignoring these interactions results in poor causal
estimates (Halloran and Struchiner, 1995; Sobel, 2006). We propose a method
that leverages the observed network structure of interactions between units to
account for treatment interference among them.

We study a setting in which a treatment has been uniformly randomized over a
set of units connected in a network, and where treatments of connected units
can influence each others’ outcomes. The development of methods for this
setting is a relatively new field in causal inference methodology, and only
few approaches for it have been proposed (e.g., van der Laan, 2014; Aronow et
al., 2017; Sussman and Airoldi, 2017).

In this paper, we propose a method that leverages matching (Rosenbaum and
Rubin, 1983) to recover direct treatment effects from experiments with
interference. Our method makes several key contributions to the study of this
setting: First, our method explicitly leverages information about the network
structure of the experimental sample to adjust for possible interference while
estimating direct treatment effects. Second, unlike other methods, matching
allows us to nonparametrically estimate treatment effects, without the need to
specify parametric models for interference or outcomes. Third, matching
produces highly interpretable results, informing analysts as to which features
of the input data were used to produce estimates. More specifically, we match
on features of graphs that are easy to interpret and visualize.

In our setting, units experience interference according to their neighborhood
graphs – the graphs defined by the units they are directly connected to. Units
with similar neighborhood graphs will experience similar interference. For
example, the educational outcome of a student randomly assigned to an extra
class depends on whether or not her friends are also assigned to that class,
and not just on how many: the specific structure of the student’s friendship
circle will influence whether or not study groups are formed, how information
is shared, how much attention the treated student will devote to the class,
and so on. All of this will impact the overall educational outcomes of
interest.

Because of this, matching units with similar neighborhood graphs together will
enable us to recover direct treatment effects even under interference. We
match units’ neighborhood graphs on counts of subgraphs within them, as graphs
with similar counts of the same unique subgraphs are naturally likely to be
similar. From there, we construct matches on individuals with similar sets of
important subgraphs; here, the set of important subgraphs is learned from a
training set. We generalize the Almost-Matching-Exactly (AME) framework (Dieng
et al., 2019; Wang et al., 2019) to match units on subgraphs in experimental
settings. We do this by constructing graph-based features that can explain
both the interference pattern in the experiment and predict the underlying
social network. We demonstrate that our method performs better than other
methods for the same problem in many settings, while generating interpretable
matches.

The paper will proceed as follows: In Section 2, we make explicit the
assumptions underpinning our framework, and outline our matching approach to
estimating direct treatment effects. In Sections 3 and 4, we evaluate the
effectiveness of our method on simulated and real-world data. Theoretical
evaluation of our approach is available in the appendix.

####  1.1 Related Work

Work on estimating causal effects under interference between units has three
broad themes. First, there has been a growing body of work on the design of
novel randomization schemes to perform causal inference under interference
(Liu and Hudgens, 2014; Sinclair et al., 2012; Duflo and Saez, 2003; Basse and
Airoldi, 2018). Some of this work makes explicit use of observed network
structure to randomly assign treatment so as to reduce interference (Ugander
et al., 2013; Toulis and Kao, 2013; Eckles et al., 2016, 2017; Jagadeesan et
al., 2019). These methodologies are inapplicable to our setting as they
require non-uniform treatment assignment, whereas in our setting we wish to
correct for interference after randomization. Second, there is work on
estimating direct treatment effects in experiments under interference, and
after conventional treatment randomization, similar to our setting. Some
existing work aims to characterize the behavior of existing estimators under
interference (Manski, 2013; Sävje et al., 2017). Other approaches lay out
methods based on randomization inference to test a variety of hypotheses under
interference and treatment randomization (Rosenbaum, 2007; Aronow, 2012; Athey
et al., 2018). Some of these approaches mix randomization inference and
outcome models (Bowers et al., 2013). For the explicit problem of recovery of
treatment effects under interference, Aronow et al. (2017) provide a general
framework to translate different assumptions about interference into inverse-
probability estimators, and Sussman and Airoldi (2017) give linearly unbiased,
minimum integrated-variance estimators under a series of assumptions about
interference. These methods either ignore explicit network structure, or
require probabilities under multiple complex sampling designs to be estimated
explicitly. Finally, there have been studies of observational inference under
network interference (van der Laan, 2014; Liu et al., 2016; Ogburn et al.,
2017; Forastiere et al., 2016). However, recovering causal estimates using
observational data when units are expected to influence each other requires a
structural model of both the nature of interference and contagion among units.

###  2 METHODOLOGY

We discuss our problem and approach in this section.

####  2.1 Problem Statement

We have a set of $n$ experimental units indexed by $i$. These units are
connected in a known graph $G=(V,E)$, where $V(G)=\\{1,\dots,n\\}$ is the set
of vertices of $G$, and $E(G)$ is the set of edges of $G$. We disallow self-
loops in our graph. We say that $H$ is a subgraph of $G$ if $V(H)\subseteq
V(G)$ and $E(H)\subseteq E(G)$. Let $t_{i}\in\\{0,1\\}$ represent the
treatment indicator for unit $i$, $\mathbf{t}$ represent the vector of
treatment indicators for the entire sample, and $\mathbf{t}_{-i}$ represent
the treatment indicators for all units except $i$. Given a treatment vector
$\mathbf{t}$ on the entire sample (i.e., for all vertices in $G$), we use
$G^{\mathbf{t}}$ to denote the labeled graph, where each vertex $i\in V(G)$
has been labeled with its treatment indicator $t_{i}$. In addition, we use
$G_{P}$ to denote a graph induced by the set of vertices $P\subseteq V(G)$ on
$G$, such that $V(G_{P})=P$ and $E(G_{P})=\\{(e_{1},e_{2})\in E(G):\;e_{1}\in
P,e_{2}\in P\\}$. We use the notation $\mathcal{N}_{i}=\\{j:(i,j)\in E(G)\\}$
to represent the neighborhood of vertex $i$. The labeled neighborhood graph of
a unit $i$, $G_{\mathcal{N}_{i}}^{\mathbf{t}}$, is defined as the graph
induced by the neighbors of $i$, and labeled according to $\mathbf{t}$. We
also define $\mathbf{t}_{\mathcal{N}_{i}}$ to be the vector of treatment
indicators corresponding to unit $i$’s neighborhood graph. A unit’s response
to the treatment is represented by its random potential outcomes
$Y_{i}(\mathbf{t})$ = $Y_{i}(t_{i},\mathbf{t}_{-i})$. Unlike other commonly
studied causal inference settings, unit $i$’s potential outcomes are now a
function of both the treatment assigned to $i$, and of all other units’
treatments. Observed treatments for unit $i$ and the whole sample are
represented by the random variables $T_{i}$ and $\mathbf{T}$ respectively. We
assume that the number of treated units is always $n^{(1)}$, i.e.,
$\sum_{i=1}^{n}T_{i}=n^{(1)}$.

A0: Ignorability of Treatment Assignment. We make the canonical assumption
that treatments are administered independently of potential outcomes, that is:
$Y_{i}(t_{i},\mathbf{t}_{-i})\mathrel{\text{\scalebox{1.07}{$\perp\mkern-10.0mu\perp$}}}\mathbf{T}$,
and $0<\Pr(T_{i}=1)<1$ for all units. In practice, we assume that treatment is
assigned uniformly at random to units, which is possible only in experimental
settings. As stated before, we do not make the canonical Stable Unit Treatment
Value Assumption (SUTVA) (Rubin, 1980), which, among other requirements,
states that units are exclusively affected by the treatment assigned to them.
We do not make this assumption because our units are connected in a network:
it could be possible for treatments to spread along the edges of the network
and to affect connected units’ outcomes. We do maintain the assumption of
comparable treatments across units, which is commonly included in SUTVA.

Our causal quantity of interest will be the Average Direct Effect (ADE), which
is defined as follows:

$\displaystyle ADE$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[Y_{i}(1,\mathbf{0})-Y_{i}(0,\mathbf{0})],$
(1)

where $\mathbf{t}_{-i}=\mathbf{0}$ represents the treatment assignment in
which no unit other than $i$ is treated. The summand represents the treatment
effect on unit $i$ when no other unit is treated, and, therefore, no
interference occurs (Halloran and Struchiner, 1995).

####  2.2 Framework

We outline the requirements of our framework for direct effect estimation
under interference. We denote interference effects on a unit $i$ with the
function $f_{i}(\mathbf{t}):\\{0,1\\}^{n}\mapsto\mathbb{R}$, a function that
maps each possible treatment allocation for the $n$ units to the amount of
interference on unit $i$. We will use several assumptions to restrict the
domain of $f$ to a much smaller set (and overload the notation $f_{i}$
accordingly). To characterize $f$, we rely on the typology of interference
assumptions introduced by Sussman and Airoldi (2017). The first three
assumptions (A0-A2) needed in our framework are common in the interference
literature (e.g., Manski, 2013; Toulis and Kao, 2013; Eckles et al., 2016;
Athey et al., 2018):

A1: Additivity of Main Effects. First, we assume that main treatment effects
are additive, i.e., that there is no interaction between units’ treatment
indicators. This allows us to write:

$\displaystyle
Y_{i}(t,\mathbf{t}_{-i})=t\tau_{i}+f_{i}(\mathbf{t}_{-i})+\epsilon_{i}$ (2)

where $\tau_{i}$ is the direct treatment effect on unit $i$, and
$\epsilon_{i}$ is some baseline effect.
A2: Neighborhood Interference. We focus on a specific form of the interference
function $f_{i}$ by assuming that the interference experienced by unit $i$
depends only on treatment of its neighbors. That is, if for two treatment
allocations $\mathbf{t},\mathbf{t}^{\prime}$ we have
$\mathbf{t}_{\mathcal{N}_{i}}=\mathbf{t}_{\mathcal{N}_{i}}^{\prime}$ then
$f_{i}(\mathbf{t})=f_{i}(\mathbf{t}^{\prime})$. To make explicit this
dependence on the neighborhood subgraph, we will write
$f_{i}(\mathbf{t}_{\mathcal{N}_{i}})\equiv f_{i}(\mathbf{t})$.
A3: Isomorphic Graph Interference We assume that, if two units $i$ and $j$
have _isomorphic labeled neighborhood graphs_ , then they receive the same
amount of interference, denoting isomorphism by $\simeq$,
$G_{\mathcal{N}_{i}}^{\mathbf{t}}\simeq
G_{\mathcal{N}_{j}}^{\mathbf{t}}\implies
f_{i}(\mathbf{t}_{\mathcal{N}_{i}})=f_{j}(\mathbf{t}_{\mathcal{N}_{j}})\equiv
f(G_{\mathcal{N}_{i}}^{\mathbf{t}})=f(G_{\mathcal{N}_{j}}^{\mathbf{t}})$.
While Assumptions A1 and A2 are standard, A3 is new. This assumption allows us
to study interference in a setting where units with similar neighborhood
subgraphs experience similar amounts of interference.

All our assumptions together induce a specific form for the potential
outcomes, namely that they depend on neighborhood structure
$G_{\mathcal{N}_{i}}^{\mathbf{t}}$, but not exactly who the neighbors are
(information contained in $\mathcal{N}_{i}$) nor treatment assignments for
those outside the neighborhood (information contained in
$\mathbf{t}_{\mathcal{N}_{i}}$). Namely:

######  Proposition 1.

Under assumptions A0-A3, potential outcomes in (2) for all units $i$ can be
written as:

$\displaystyle Y_{i}(t,\mathbf{t}_{-i})$
$\displaystyle=t\tau_{i}+f(G_{\mathcal{N}_{i}}^{\mathbf{t}})+\epsilon_{i},$
(3)

where $\tau_{i}$ is the direct treatment effect on unit $i$, and
$\epsilon_{i}$ is some baseline response.

In addition, suppose that baseline responses for all units are equal to each
other in expectation, i.e., for all $i$, $\mathbb{E}[\epsilon_{i}]=\alpha$.
Then under assumptions A0-A3, for neighborhood graph structures $g_{i}$ of
unit $i$ and treatment vectors $\mathbf{t}$, the ADE is identified as:

$\displaystyle
ADE=\frac{1}{n^{(1)}}\sum_{i=1}^{n}\mathbb{E}\bigl{[}T_{i}\times$
$\displaystyle\bigl{(}\mathbb{E}[Y_{i}|G_{\mathcal{N}_{i}}^{\mathbf{T}}\simeq
g_{i}^{\mathbf{t}},T_{i}=1]$
$\displaystyle-\mathbb{E}[Y_{i}|G_{\mathcal{N}}^{\mathbf{T}}\simeq
g_{i}^{\mathbf{t}},T_{i}=0]\bigr{)}\bigr{]},$

where $G_{\mathcal{N}_{i}}^{\mathbf{T}}$ is the neighborhood graph of $i$
labelled according to the treatment assignment $\mathbf{T}$.

The proposition (whose proof is in the appendix) states that the interference
received by a unit is a function of each unit’s neighborhood graph. Further,
the outcomes can be decomposed additively into this function and the direct
treatment effect on $i$. The proposition implies that the ADE is identified by
matching each treated unit to one or more control units with an isomorphic
neighborhood graph, and computing the direct effect on the treated using these
matches. This effect is, in expectation over individual treatment assignments,
equal to the ADE.

####  2.3 Subgraph Selection via Almost-Matching-Exactly

Given Proposition 1 and the framework established in the previous section, we
would ideally like to match treated and control units that have isomorphic
neighborhood graphs. This would allow us to better estimate the ADE without
suffering interference bias: for a treated unit $i$, if a control unit $j$ can
be found such that $G_{\mathcal{N}_{i}}^{\mathbf{t}}\simeq
G_{\mathcal{N}_{j}}^{\mathbf{t}}$, then $j$’s outcome will be identical in
expectation to $i$’s counterfactual outcome and can be used as a proxy.
Unfortunately, the number of non-isomorphic (canonically unique) graphs with a
given number of nodes and edges grows incredibly quickly (Harary, 1994) and
finding such matches is infeasible for large graphs. We therefore resort to
counting all subgraphs that appear in a unit’s neighborhood graph and matching
units based on the counts of those subgraphs. However, instead of exactly
matching on the counts of those subgraphs, we match treated and control units
if they have _similar_ counts, since matching exactly on all subgraph counts
implies isomorphic neighborhoods and is also infeasible. Further, absolutely
exact matches may not exist in real networks.

Constructing inexact matches, in turn, requires a measure of relative graph
importance. In Figure 1, for example, there are two control units that the
treated unit may be matched to; if triangles contribute more to the
interference function, it should be matched to the right; otherwise, if degree
and/or two-stars are more important, it should be matched to the left. Of
course, these relative importance measures might depend on the problem and we
would like to learn them.

Figure 1: Inexact matching presupposes an ordering of feature importance;
should the the treated ego (black) be matched to a control whose neighborhood
graph has the same number of units (left), or same number of triangles
(right)?

It might be tempting to match directly on $f$, as that would lead to unbiased
inference. However, we abstain from doing so for two reasons. Firstly, in
practice, the true interference is unknown and we could only match on
estimated values of $f$; this suffers from all the problems that afflict
matching on estimated propensity scores without appropriate adjustments
(Abadie and Imbens, 2016) or parametric approximations (Rubin and Thomas,
1996). Such corrections or approximations do not currently exist for estimated
interference functions and their development is an active area of research.
Secondly, interpretability is a key component of our framework that would be
lost matching on $f$-values; these values are scalar summaries of interference
that depends on entire graphs. Estimating $f$ well would also likely require
complex and uninterpretable nonparametric methods. In Section K of the
appendix, we empirically compare matching units on $f$-values to our subgraph
matching method via simulation. The loss of interpretability associated with
matching on $f$ does not yield substantial gains in performance, even when
using _true_ values of $f$ for matching, which is impossible in practice.

Almost-Matching-Exactly (AME) (Wang et al., 2019; Dieng et al., 2019; Awan et
al., 2019) provides a framework for the above problem that is explicitly
geared towards building interpretable, high-quality matches on discrete
covariates, which in our setting are the counts of the treated subgraphs in
the neighborhood. AME performs inexact matching while learning importance
weights for each covariate from a training set, prioritizing matches on more
important covariates. In this way, it neatly addresses the challenge of
inexact matching by learning a metric specific to discrete covariates (namely,
a weighted Hamming distance). Formally, AME matches units so as to optimize a
flexible measure of match quality. For each treated unit $i$, solving the AME
problem is equivalent to finding:

$\displaystyle\boldsymbol{\theta}^{i^{*}}$
$\displaystyle\in\operatorname*{arg\,max}_{\boldsymbol{\theta}\in\\{0,1\\}^{p}}\boldsymbol{\theta}^{T}\mathbf{w}$
(4) $\displaystyle\text{such that }\exists j:t_{j}=0\text{ and
}\mathbf{x}_{j}\circ\boldsymbol{\theta}=\mathbf{x}_{i}\circ\boldsymbol{\theta}$

where $\circ$ denotes the Hadamard product, $\mathbf{w}$ is a vector of
weights and $\boldsymbol{x}_{i},\boldsymbol{x}_{j}$ are vectors of binary
covariates for units $i$ and $j$ that we might like to match on. In our
network interference setting, these are vectors of subgraph counts. The vector
$\mathbf{w}$ denotes the importance of each subgraph in causing interference.
We will leverage both information on outcomes and networks to construct an
estimate for it.

We start by enumerating (up to isomorphism) all the $p$ subgraphs
$g_{1},\dots,g_{p}$ that appear in any of the
$G_{\mathcal{N}_{i}}^{\mathbf{t}},i\in 1,\dots,n$. The covariates for unit $i$
are then given by
$S(G_{\mathcal{N}_{i}}^{\mathbf{t}})=(S_{1}(G_{\mathcal{N}_{i}}^{\mathbf{t}}),\dots,S_{p}(G_{\mathcal{N}_{i}}^{\mathbf{t}}))$
where $S_{k}(G_{\mathcal{N}_{i}}^{\mathbf{t}})$ denotes the number of times
subgraph $g_{k}$ appears in the subgraphs of
$G_{\mathcal{N}_{i}}^{\mathbf{t}}$. These counts are then converted into
binary indicators that are one if the count of subgraph $g_{k}$ in each unit’s
neighborhood is exactly $x$, for all $x$ observed in the data. Thus, units
will be matched exactly if they have identical subgraph counts. We then
approximately solve the problem in Equation (4) to find the optimally
important set of subgraphs upon which to exactly match each treated unit, such
that there is at least one control unit that matches exactly with the treated
unit on the chosen subgraph counts. The key idea behind this approach is that
we want to match units exactly on subgraph counts that contribute
significantly to the interference function, trading off exactly-matching on
these important subgraphs with potential mismatches on subgraphs that
contribute less to interference.

In practice, our implementation enumerates all subgraphs in each unit’s
neighborhood and stores the count of each pattern – this is computationally
challenging. There is a growing body of work on efficient counting algorithms
for pre-specified small patterns (up to 4-5 nodes) but there is little
research on fast methods to both enumerate and count all motifs in a graph
(e.g., Pinar et al., 2017; Marcus and Shavitt, 2010; Hu et al., 2013).
Empirically, we see that this enumeration takes less than 30 seconds for 50
units.

###### The FLAME Algorithm for AME.

The Fast Large Almost Matching Exactly (FLAME) algorithm (Wang et al., 2019)
approximates the solution to the AME problem. The procedure starts by exactly
matching all possible units on all covariates. It then drops one covariate at
a time, choosing the drop maximizing the match quality $\mathtt{MQ}$ at that
iteration, defined:

$\displaystyle\mathtt{MQ}=C\cdot\mathtt{BF}-\widehat{\mathtt{PE}}_{Y}.$ (5)

The match quality is the sum of a balancing factor $\mathtt{BF}$ and a
predictive error $\widehat{\mathtt{PE}}_{Y}$, with relative weights determined
by the hyper-parameter $C$. The balancing factor is defined as the proportion
of treated units plus the proportion of control units matched at that
iteration. Introducing the balancing factor into the objective has the
advantage of encouraging more units to be matched, thereby minimizing variance
of estimators (see Wang et al., 2019). In our setting, the second component of
the match quality, predictive error, takes the form:

$\displaystyle\widehat{\mathtt{PE}}_{Y}=\operatorname*{arg\,min}_{h\in\mathcal{F}_{1}}\sum_{i}^{n}(Y_{i}-h(S(G_{\mathcal{N}_{i}}^{\mathbf{t}})\circ\boldsymbol{\theta},T_{i}))^{2}$
(6)

for some class of functions $\mathcal{F}_{1}$. It is computed using a holdout
training set and discourages dropping covariates that are useful for
predicting the outcome. In this way, FLAME strikes a balance between matching
many units and ensuring these matches are of high-quality. By using a holdout
set to determine how useful a set of variables is for out-of-sample
prediction, FLAME learns a measure of covariate importance via a weighted
Hamming distance. Specifically, it learns a vector of importance weights
$\mathbf{w}$ for the different subgraph counts that minimizes
$\mathbf{w}^{T}\mathbb{I}[S(G_{\mathcal{N}_{i}}^{\mathbf{t}})\neq
S(G_{\mathcal{N}_{j}}^{\mathbf{t}})]$, where
$\mathbb{I}[S(G_{\mathcal{N}_{i}}^{\mathbf{t}})\neq
S(G_{\mathcal{N}_{j}}^{\mathbf{t}})]$ is a vector whose $k^{\textrm{th}}$
entry is 0 if the labeled neighborhood graphs of $i$ and $j$ have the same
count of subgraph $k$, and 1 otherwise.

To this match quality term, we add a network fit term to give subgraphs more
weight that are highly predictive of overall network structure. We fit a
logistic regression model in which the edges $(i,j)$ between units $i,j$ are
independent given $\mathcal{N}_{i},\mathcal{N}_{j}$, and dependent on the
subgraph counts of units $i$ and $j$:

$(i,j)\stackrel{{\scriptstyle
iid}}{{\sim}}\textrm{Bern}(\textrm{logit}(\beta_{1}^{T}S(G_{\mathcal{N}_{i}}^{\mathbf{t}})+\beta_{2}^{T}S(G_{\mathcal{N}_{j}}^{\mathbf{t}})))$

To the match quality in the original formulation, we then add
$\widehat{\mathtt{PE}}_{G}$, defined to be the _AIC_ (Akaike, 1974) of this
fitted model, weighted by a hyperparameter $D$. Therefore, at each iteration:

$\displaystyle\mathtt{MQ}=C\cdot\mathtt{BF}-\widehat{\mathtt{PE}}_{Y}+D\cdot\widehat{\mathtt{PE}}_{G}.$

Thus, we penalize not only subgraph drops that impede predictive performance
or making matches, but also those that make the observed network unlikely.
$\widehat{\mathtt{PE}}_{G}$ represents the empirical prediction error of the
chosen set of statistics for the observed graph: if
$\widehat{\mathtt{PE}}_{G}$ is low, then the chosen subgraphs do a good job of
predicting the observed graph. This error is also evaluated at a minimum over
another class of prediction functions, $\mathcal{F}_{2}$. This term in the AME
objective is justified by Assumption A3: units that have isomorphic labeled
neighborhood graphs should experience the same amount of interference, and
subgraph counts should be predictive of neighborhood graph structure.

Our approach to estimating the ADE is therefore as follows. (1) For each unit
$i$, count and label all types of subgraphs in
$G_{\mathcal{N}_{i}}^{\mathbf{t}}$. (2) Run FLAME, encouraging large numbers
of matches on subgraph counts, while using the covariates that are most
important for predicting the outcome and the network. (3) Estimate ADE as
$\widehat{ADE}$, by computing the difference in means for each matched group
and then averaging across matched groups, weighted by their size. Since our
approach is based on FLAME, we call it _FLAME-Networks_.

Extensions. FLAME-Networks immediately extends to handling unit-level
covariate information for baseline adjustments; we simply concatenate subgraph
information and covariate information in the same dataset and then make
almost-exact matches. FLAME-Networks will automatically learn the weights of
both subgraphs and baseline covariates to make matches that take both into
account. Another straightforward extension considers interference not just in
the immediate neighborhood of each unit, but up to an arbitrary number of hops
away. To extend FLAME-Networks in this way, it is sufficient to enumerate
subgraphs in the induced neighborhood graph of each unit $i$ where the
vertices considered are those with a path to $i$ that is at most $k$ steps
long. Given these counts, our method proceeds exactly as before.

###  3 EXPERIMENTS

We begin by evaluating the performance of our estimator in a variety of
simulated settings, in which we vary the form of the interference function. We
find that our approach performs well in many distinct settings. In Section M
of the appendix, we also assess the quality of the matches constructed.

We simulate graphs from an Erdős-Rènyi model:
$G\sim\textrm{Erd\H{o}s-R\\`{e}nyi}(n,q)$, by which every possible edge
between the $n$ units is created independently, with probability $q$. In
Sections H and I of the appendix, we also perform experiments on cluster-
randomized and real-world networks. Treatments for the whole sample are
generated with $\Pr(\mathbf{T}=\mathbf{t})=\binom{n}{n^{(1)}}^{-1}$, where
$n^{(1)}$ is the number of treated units. Outcomes are generated according to
$Y_{i}(t,\mathbf{t}_{-i})=\mathbf{t}\tau_{i}+f(G_{\mathcal{N}_{i}}^{\mathbf{t}})+\epsilon_{i}$,
where $\epsilon\sim N(\mathbf{0},I_{n})$ represents a baseline outcome;
$\tau_{i}\sim N(\mathbf{5},I_{n})$ represents a direct treatment effect, and
$f$ is the interference function. In Section J of the appendix, we consider a
setting in which the errors are heteroscedastic. For the interference
function, we use additive combinations of the subgraph components in Table 1
and define $m_{ip},p=1,\dots,7$ to be the counts of feature $p$ in
$G_{\mathcal{N}_{i}\cup\\{i\\}}^{\mathbf{t}}$. Lastly, the counts of each
component are normalized to have mean 0 and standard deviation 1.

$d_{i}$: |  The treated degree of unit $i$
---|---
$\Delta_{i}$: |  The number of triangles in $G_{\mathcal{N}_{i}\cup\\{i\\}}^{\mathbf{t}}$ with at least one treated unit.
$\bigstar_{i}^{k}$: |  The number of $k$-stars in $G_{\mathcal{N}_{i}\cup\\{i\\}}^{\mathbf{t}}$ with at least one treated unit.
$\dagger_{i}^{k}$: |  The number of units in $G_{\mathcal{N}_{i}\cup\\{i\\}}^{\mathbf{t}}$ with degree $\geq k$ and at least one treated unit among their neighbors.
$B_{i}$: |  The vertex betweenness of unit $i$.
$C_{i}$: |  The closeness centrality of unit $i$.
Table 1: Interference components used in experiments; see the appendix for
more details.

We compare our approach with different methods to estimate the ADE under
interference:

Naïve. The simple difference in means between treatment and control groups
assuming no interference.
All Eigenvectors. Eigenvectors for the entire adjacency matrix are computed
with every treated unit matched to the control unit minimizing the Mahalanobis
distance between the eigenvectors, weighing the $k$’th eigenvector by $1/k$.
The idea behind this estimator is that the eigendecomposition of the adjacency
matrix encodes important information about the network and how interference
might spread within it.
First Eigenvector. Same as All Eigenvectors except units are matched only on
their values of the largest-eigenvalue eigenvector.
Stratified Naïve. The stratified naïve estimator as discussed by Sussman and
Airoldi (2017). A weighted difference-in-means estimator where units are
divided into strata defined by their treated degree (number of treated
vertices they are connected to), and assigned weight equal to the number of
units within the stratum in the final difference of weighted averages between
treated and control groups.
SANIA MIVLUE. The minimum integrated variance, linear unbiased estimator under
assumptions of symmetrically received interference and additivity of main
effects, when the priors on the baseline outcome and direct treatment effect
have no correlation between units; proposed by Sussman and Airoldi (2017).
FLAME-Networks. Our proposed method. In all simulations, the two components of
the PE function are weighted equally, and a ridge regression is used to
compute outcome prediction error.

####  3.1 Experiment 1: Additive Interference

First we study a setting in which interference is an additive function of the
components in Table 1. Outcomes in this experiment have the form:
$Y_{i}=\gamma_{1}d_{i}+\gamma_{2}\Delta_{i}+\gamma_{3}\bigstar^{2}_{i}+\gamma_{4}\bigstar^{4}_{i}+\gamma_{5}\dagger_{i}^{3}+\gamma_{6}B_{i}+\gamma_{7}C_{i}+\epsilon_{i}$,
with $\epsilon_{i}\sim N(0,1)$. We simulate 50 datasets for each setting, in
which the units are in an $ER(50,0.05)$ graph. Table 2 shows values for the
$\gamma_{i}$ in each of our experimental settings.

Feature | $d_{i}$ | $\Delta_{i}$ | $\bigstar^{2}_{i}$ | $\bigstar^{4}_{i}$ | $\dagger^{3}$ | $B_{i}$ | $C_{i}$
---|---|---|---|---|---|---|---
Weight | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$ | $\gamma_{4}$ | $\gamma_{5}$ | $\gamma_{6}$ | $\gamma_{7}$
Setting 1 | 0 | 10 | 0 | 0 | 0 | 0 | 0
Setting 2 | 10 | 10 | 0 | 0 | 0 | 0 | 0
Setting 3 | 0 | 10 | 1 | 1 | 1 | 1 | -1
Setting 4 | 5 | 1 | 10 | 1 | 1 | 1 | -1
Table 2: Settings for Experiment 1.

Results for Experiment 1 are reported in Figure 2. FLAME-Networks outperforms
all other methods both in terms of average error, and standard deviation over
the simulations. This is likely because FLAME-Networks learns weights for the
subgraphs that are proportionate to those we use at each setting, and matches
units on subgraphs with larger weights. When the interference function is
multiplicative instead of additive, FLAME-Networks performs similarly; results
are in the appendix.

Figure 2: Results from Experiment 1. Each violin plot represents the
distribution over simulations of absolute estimation error. The panels are
numbered according to the parameter settings of the simulations. Violin plots
are blue if the method had mean error lower than or equal to FLAME-Networks’
and red otherwise. The black line inside each violin is the median error. The
dashed line is FLAME-Networks’ mean error.

####  3.2 Experiment 2: Covariate Adjustment

A strength of FLAME-Networks is its ability to natively account for covariate
information. We analyze a setting in which baseline effects are dependent on
an additional discrete-valued covariate, $x$, that is observed alongside the
network. Outcomes take the form
$Y_{i}=t\tau_{i}+f(G_{\mathcal{N}_{i}}^{\mathbf{t}})+\beta
x_{i}+\epsilon_{i}$, where $x_{i}$ is chosen uniformly at random from
$\\{1,2,3\\}$ for each unit, and $\beta$ is fixed at 15. This means that our
sample is divided into 3 strata defined by the covariate values. We ran FLAME-
Networks with a dataset consisting of subgraph counts, plus observed
covariates for each unit. For comparison with the other methods, we first
regress $Y$ on $x$, and then use the residuals of that regression as outcomes
for the other methods. This way, the initial regression will account for the
baseline effects of $x$, and the residuals contain only potential
interference. The interference function takes the form
$f(G_{\mathcal{N}_{i}}^{\mathbf{t}})=d_{i}+\Delta_{i}+B_{i}$, which is what we
are trying to learn with the other methods. We simulate the sample network
from $ER(50,0.05)$.

Results are displayed in Table 3. FLAME-Networks performs, on average, better
than all the other methods. Results in Section L of the supplement show that
when $\beta$ is increased, none of the methods suffer in performance. While
regression adjustment prior to estimation seems to have a positive impact on
the performance of other methods in the presence of additional covariates,
FLAME-Networks performs best. This is because FLAME-Networks is built to
easily handle the inclusion of covariates in its estimation procedure.

Method | Median | 25th q | 75th q
---|---|---|---
FLAME-Networks | 0.39 | 0.21 | 0.59
First Eigenvector | 0.47 | 0.40 | 0.83
All Eigenvectors | 0.55 | 0.29 | 0.79
Naive | 0.53 | 0.36 | 0.92
SANIA | 1.93 | 1.75 | 2.25
Stratified | 4.49 | 4.45 | 4.53
Table 3: Results from Experiment 2 with $\beta=5$. Median and 25th and 75th
percentile of absolute error over 40 simulated datasets.

####  3.3 Experiment 3: Misspecified Interference

We now study the robustness of our estimator in a setting in which one of our
key assumptions – A3 – is violated. Specifically, we now allow for treated and
control units to receive different amounts of interference, even if they have
the same labelled neighborhood graphs. We do this by temporarily eliminating
all control-control edges in the network and then counting the features in
Table 1 used to assign interference. That is, consider a unit $i$ with a
single, untreated neighbor $j$. In our new setting, if degree is a feature of
the interference function, then $i$ being treated implies $i$ receives
interference from $j$. But if $i$ is untreated, then $i$ would receive no
interference from $j$, because its neighbor is also untreated. This crucially
implies that FLAME-Networks will be matching–and estimating the ADE
from–individuals that do not necessarily receive similar amounts of
interference.

In this setting, we generate interference according to:
$f_{i}=(5-\gamma)d_{i}+\gamma\Delta_{i}$ for $\gamma\in[0,5]$ and assess the
performance of FLAME-Networks against that of the SANIA and stratified
estimators. Results are shown in Figure 3. We see that, when degree is the
only component with weight in the interference function, FLAME-Networks
performs better than the stratified estimator, but worse than the SANIA
estimator, which leverages aspects of the graph related to degree. However,
our performance improves as $\gamma$ increases and the true interference
depends more on triangle counts since the triangle counts available to FLAME-
Networks represent the actual interference pattern more frequently than the
degree counts did. Thus, we see that although violation of our method’s
assumptions harms its performance, it still manages at times to outperform
estimators that rely too heavily on degree.

Figure 3: Results from Experiment 3. These simulations were run on an ER(75,
0.07) graph. The bands around the lines represent 25th and 75th quantiles of
the 50 simulations, for each value of $\gamma$.

###  4 APPLICATION

In this section, we demonstrate the practical utility of our method. We use
data collected by Banerjee et al. (2013) on social networks for 75 villages in
Karnataka, India. They are a median distance of 46 km from one another other,
motivating the assumption that network interference is experienced solely
between individuals from the same village. For each village, we study the
effect of 1. lack of education on election participation; 2. lack of education
on Self-Help-Group (SHG) participation; and 3. being male on SHG
participation. We proxy election participation by ownership of an election
card. We compare our estimates – which account for network interference – to
naive estimates – which assume no network interference. Data pre-processing is
summarized in the appendix.

For ADE estimates, we assume the treatment is randomly assigned. We find that
lack of education is associated with higher SHG participation, and that males
are less likely to participate in SHGs than females (see Figure 4). These
results make sense in the context of developing countries where SHGs are
mainly utilized by females in low-income families, which generally have lower
education levels. We observe that education does not impact election
participation; in developing countries, an individual’s decision to
participate in an election may be driven by factors such caste, religion,
influence of local leaders and closeness of race (Shachar and Nalebuff, 1999;
Gleason, 2001). FLAME-Networks matches units in each village by subgraph
counts and covariate values to estimate the ADE. Looking at the matched
groups, we discover that subgraphs such as 2-stars and triangles were
important for matching, implying that second-order connections could be
affecting interference in this setting. Further details of the matched groups
are in Section F.

Figure 4 plots naive and FLAME-Networks ADE estimates. We find a significant
difference between our estimates and the naive estimates when estimating the
effect of being male on participation in SHGs. The naive estimator
overestimates the treatment effect, which is possible when ignoring network
interference. It is plausible, in this setting, that interference would
heighten these effects for all outcomes. This is because individuals from
similar social backgrounds or gender tend to interact more together, and,
therefore, are more likely to influence each other’s participation decision,
both in elections and in SHGs.

Figure 4: Naive and FLAME-Networks ADE estimates and their difference. Red,
blue, and green respectively correspond to (treatment, outcome) pairs: (no
education, election participation), (no education, SHG participation), and
(gender, SHG participation).

###  5 DISCUSSION

Conventional estimators for treatment effects in randomized experiments will
be biased when there is interference between units. We have introduced FLAME-
Networks – a method to recover direct treatment effects in such settings. Our
method is based on matching units with similar neighborhood graphs in an
almost-exact way, thus producing interpretable, high-quality results. We have
shown that FLAME-Networks performs better than existing methods both on
simulated data, and we have used real-world data to show how it can be applied
in a real setting. Our method extends easily to settings with additional
covariate information for units and to taking into account larger
neighborhoods for interference. In future work, our method can be extended to
learning a variety of types of distance metrics between graphs.

##### Acknowledgements

This work was supported in part by NIH award 1R01EB025021-01, NSF awards
IIS-1552538 and IIS1703431, a DARPA award under the L2M program, and a Duke
University Energy Initiative ERSF grant.

### References

  * Abadie and Imbens (2016) Alberto Abadie and Guido Imbens.  Matching on the estimated propensity score.  _Econometrica_ , pages 781–807, 2016.
  * Akaike (1974) H. Akaike.  A new look at the statistical model identification.  _IEEE Transactions on Automatic Control_ , 19(6):716–723, December 1974.  doi: 10.1109/TAC.1974.1100705.
  * Aronow (2012) Peter M Aronow.  A general method for detecting interference between units in randomized experiments.  _Sociological Methods & Research_, 41(1):3–16, 2012.
  * Aronow et al. (2017) Peter M Aronow, Cyrus Samii, et al.  Estimating average causal effects under general interference, with application to a social network experiment.  _The Annals of Applied Statistics_ , 11(4):1912–1947, 2017.
  * Athey et al. (2018) Susan Athey, Dean Eckles, and Guido W Imbens.  Exact p-values for network interference.  _Journal of the American Statistical Association_ , 113(521):230–240, 2018.
  * Awan et al. (2019) M Awan, Yameng Liu, Marco Morucci, Sudeepa Roy, Cynthia Rudin, and Alexander Volfovsky.  Interpretable almost matching exactly with instrumental variables.  In _Conference on Uncertainty in Artificial Intelligence (UAI)_ , 2019\.
  * Banerjee et al. (2013) Abhijit Banerjee, Arun G Chandrasekhar, Esther Duflo, and Matthew O Jackson.  The diffusion of microfinance.  _Science_ , 341(6144):1236498, 2013.
  * Basse and Airoldi (2018) Guillaume W Basse and Edoardo M Airoldi.  Model-assisted design of experiments in the presence of network-correlated outcomes.  _Biometrika_ , 105(4):849–858, 2018.
  * Bowers et al. (2013) Jake Bowers, Mark M Fredrickson, and Costas Panagopoulos.  Reasoning about interference between units: A general framework.  _Political Analysis_ , 21(1):97–124, 2013.
  * Cox (1958) David Roxbee Cox.  _Planning of Experiments_.  Wiley, 1958.
  * Dieng et al. (2019) Awa Dieng, Yameng Liu, Sudeepa Roy, Cynthia Rudin, and Alexander Volfovsky.  Interpretable almost-exact matching for causal inference.  In _Proceedings of Artificial Intelligence and Statistics (AISTATS)_ , pages 2445–2453, 2019.
  * Duflo and Saez (2003) Esther Duflo and Emmanuel Saez.  The role of information and social interactions in retirement plan decisions: Evidence from a randomized experiment.  _The Quarterly Journal of Economics_ , 118(3):815–842, 2003.
  * Eckles et al. (2016) Dean Eckles, René F Kizilcec, and Eytan Bakshy.  Estimating peer effects in networks with peer encouragement designs.  _Proceedings of the National Academy of Sciences_ , 113(27):7316–7322, 2016.
  * Eckles et al. (2017) Dean Eckles, Brian Karrer, and Johan Ugander.  Design and analysis of experiments in networks: Reducing bias from interference.  _Journal of Causal Inference_ , 5(1), 2017.
  * Forastiere et al. (2016) Laura Forastiere, Edoardo M Airoldi, and Fabrizia Mealli.  Identification and estimation of treatment and interference effects in observational studies on networks.  _arXiv preprint arXiv:1609.06245_ , 2016.
  * Gleason (2001) Suzanne Gleason.  Female political participation and health in india.  _The ANNALS of the American Academy of Political and Social Science_ , 573(1):105–126, 2001.
  * Halloran and Struchiner (1995) M Elizabeth Halloran and Claudio J Struchiner.  Causal inference in infectious diseases.  _Epidemiology (Cambridge, Mass.)_ , 6(2):142–151, 1995.
  * Harary (1994) Frank Harary.  _Graph Theory_.  Addison-Wesley, 1994.
  * Harris (2013) Kathleen Harris.  The add health study: design and accomplishments.  Technical report, Carolina Population Center: University of North Carolina at Chapel Hill, 2013.
  * Hu et al. (2013) Xiaocheng Hu, Yufei Tao, and Chin-Wan Chung.  Massive graph triangulation.  In _Proceedings of the 2013 ACM SIGMOD international conference on Management of data_ , pages 325–336, 2013.
  * Jagadeesan et al. (2019) Ravi Jagadeesan, Natesh Pillai, and Alexander Volfovsky.  Designs for estimating the treatment effect in networks with interference.  _Annals of Statistics_ , 2019.
  * Liu and Hudgens (2014) Lan Liu and Michael G Hudgens.  Large sample randomization inference of causal effects in the presence of interference.  _Journal of the American Statistical Association_ , 109(505):288–301, 2014.
  * Liu et al. (2016) Lan Liu, Michael G Hudgens, and Sylvia Becker-Dreps.  On inverse probability-weighted estimators in the presence of interference.  _Biometrika_ , 103(4):829–842, 2016.
  * Manski (2013) Charles F Manski.  Identification of treatment response with social interactions.  _The Econometrics Journal_ , 16(1):S1–S23, 2013\.
  * Marcus and Shavitt (2010) Dror Marcus and Yuval Shavitt.  Efficient counting of network motifs.  In _2010 IEEE 30th International Conference on Distributed Computing Systems Workshops_ , pages 92–98. IEEE, 2010.
  * Ogburn et al. (2017) Elizabeth L Ogburn, Tyler J VanderWeele, et al.  Vaccines, contagion, and social networks.  _The Annals of Applied Statistics_ , 11(2):919–948, 2017.
  * Pinar et al. (2017) Ali Pinar, C Seshadhri, and Vaidyanathan Vishal.  Escape: Efficiently counting all 5-vertex subgraphs.  In _Proceedings of the 26th International Conference on World Wide Web_ , pages 1431–1440, 2017.
  * Rosenbaum (2007) Paul R Rosenbaum.  Interference between units in randomized experiments.  _Journal of the American Statistical Association_ , 102(477):191–200, 2007.
  * Rosenbaum and Rubin (1983) Paul R Rosenbaum and Donald B Rubin.  The central role of the propensity score in observational studies for causal effects.  _Biometrika_ , 70(1):41–55, 1983.
  * Rubin and Thomas (1996) Donald Rubin and Neal Thomas.  Matching using estimated propensity scores: relating theory to practice.  _Biometrics_ , 52(1):249–264, 1996.
  * Rubin (1980) Donald B Rubin.  Randomization analysis of experimental data: The fisher randomization test comment.  _Journal of the American Statistical Association_ , 75(371):591–593, 1980.
  * Sävje et al. (2017) Fredrik Sävje, Peter M Aronow, and Michael G Hudgens.  Average treatment effects in the presence of unknown interference.  _arXiv preprint arXiv:1711.06399_ , 2017.
  * Shachar and Nalebuff (1999) Ron Shachar and Barry Nalebuff.  Follow the leader: Theory and evidence on political participation.  _American Economic Review_ , 89(3):525–547, 1999\.
  * Sinclair et al. (2012) Betsy Sinclair, Margaret McConnell, and Donald P Green.  Detecting spillover effects: Design and analysis of multilevel experiments.  _American Journal of Political Science_ , 56(4):1055–1069, 2012.
  * Sobel (2006) Michael E Sobel.  What do randomized studies of housing mobility demonstrate? causal inference in the face of interference.  _Journal of the American Statistical Association_ , 101(476):1398–1407, 2006.
  * Sussman and Airoldi (2017) Daniel L Sussman and Edoardo M Airoldi.  Elements of estimation theory for causal effects in the presence of network interference.  _arXiv preprint arXiv:1702.03578_ , 2017.
  * Toulis and Kao (2013) Panos Toulis and Edward Kao.  Estimation of causal peer influence effects.  In _International Conference on Machine Learning (ICML)_ , pages 1489–1497, 2013.
  * Ugander et al. (2013) Johan Ugander, Brian Karrer, Lars Backstrom, and Jon Kleinberg.  Graph cluster randomization: Network exposure to multiple universes.  In _Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 329–337. ACM, 2013.
  * van der Laan (2014) Mark J van der Laan.  Causal inference for a population of causally connected units.  _Journal of Causal Inference_ , 2(1):13–74, 2014\.
  * Wang et al. (2019) Tianyu Wang, Marco Morucci, M Awan, Yameng Liu, Sudeepa Roy, Cynthia Rudin, and Alexander Volfovsky.  FLAME: A fast large-scale almost matching exactly approach to causal inference.  _arXiv preprint arXiv:1707.06315_ , 2019.

## Appendix

###  A Proofs from the Paper

####  A.1 Proof of Proposition 1

###### Proof.

We start by writing potential outcomes for an arbitrary unit $i$ from (2) as
$Y_{i}(t,\mathbf{t}_{-i})=f_{i}(t,\mathbf{t}_{-i})+\epsilon_{i}$, where
$\epsilon_{i}$ is some pre-treatment value of the potential outcome, and
$f_{i}(t,\mathbf{t}_{-i})$ is the treatment response function, dependent both
on unit $i$’s treatment and on everyone else’s. Using A0-A3 we can write (2)
as:

$\displaystyle Y_{i}(t,\mathbf{t}_{-i})$
$\displaystyle=t\tau_{i}+f_{i}(\mathbf{t}_{-i})+\epsilon_{i}$ (By A1)
$\displaystyle=t\tau_{i}+f_{i}(\mathbf{t}_{\mathcal{N}_{i}})+\epsilon_{i}$ (By
A2) $\displaystyle=t\tau_{i}+f(G_{\mathcal{N}_{i}}^{\mathbf{t}})+\epsilon_{i}$
$\displaystyle\mbox{(By A3)}.$ (7)

This proves Equation (3). The first line agrees with the additivity
assumption, which permits a contribution $f_{i}(\mathbf{t}_{-i})$ to
$Y_{i}(t,\mathbf{t}_{-i})$ arising from anywhere within the rest of the graph.
The second line states that this contribution is limited to $i$’s
neighborhood. The third line states that the contribution depends only the
neighborhood structure and not anything else about the neighbors.

To prove identification for the ADE, we must show that the individual
treatment effect is identified for a treated unit $i$. We need to show that:
$\mathbb{E}[Y_{i}(1,\mathbf{0})-Y_{i}(0,\mathbf{0})]=\mathbb{E}[Y_{i}|T_{i}=1,G_{\mathcal{N}_{i}}^{\mathbf{T}}=g_{i}^{\mathbf{t}}]-\mathbb{E}[Y_{i}|T_{i}=0,G_{\mathcal{N}_{i}}^{\mathbf{T}}=g_{i}^{\mathbf{t}}]$.
We have first:

$\displaystyle\mathbb{E}[Y_{i}(1,\mathbf{0})-Y_{i}(0,\mathbf{0})]$
$\displaystyle=\mathbb{E}[t\tau_{i}+f(G_{\mathcal{N}_{i}}^{\mathbf{0}})+\epsilon_{i}-f(G_{\mathcal{N}_{i}}^{\mathbf{0}})-\epsilon_{i}]$
$\displaystyle=\mathbb{E}[\epsilon_{i}+\tau_{i}-\epsilon_{i}]=\tau_{i}.$ (8)

Second, we have:

$\displaystyle\mathbb{E}[Y_{i}|G_{\mathcal{N}_{i}}^{\mathbf{T}}\simeq
g_{i}^{\mathbf{t}},T_{i}=1]-\mathbb{E}[Y_{i}|G_{\mathcal{N}_{i}}^{\mathbf{T}}\simeq
g_{i}^{\mathbf{t}},T_{i}=0]$ $\displaystyle=$
$\displaystyle\mathbb{E}[Y_{i}(1,\mathbf{T}_{-i})T_{i}+$
$\displaystyle\quad+Y_{i}(0,\mathbf{T}_{-i})(1-T_{i})|G_{\mathcal{N}_{i}}^{\mathbf{T}}\simeq
g_{i}^{\mathbf{t}},T_{i}=1]$ $\displaystyle-$
$\displaystyle\mathbb{E}[Y_{i}(1,\mathbf{T}_{-i})T_{i}$
$\displaystyle\quad+Y_{i}(0,\mathbf{T}_{-i})(1-T_{i})|G_{\mathcal{N}_{i}}^{\mathbf{T}}\simeq
g_{i}^{\mathbf{t}},T_{i}=0]$ $\displaystyle=$
$\displaystyle\mathbb{E}[\tau_{i}+f(g_{i}^{\mathbf{t}})+\epsilon_{i}|G_{\mathcal{N}_{i}}^{\mathbf{T}}\simeq
g_{i}^{\mathbf{t}},T_{i}=1]$
$\displaystyle-\mathbb{E}[f(g_{i}^{\mathbf{t}})+\epsilon_{i}|G_{\mathcal{N}_{j}}^{\mathbf{T}}=g_{i}^{\mathbf{t}},T_{i}=0]$
$\displaystyle=$
$\displaystyle\alpha+f(g_{i}^{\mathbf{t}})+\tau_{i}-\alpha-f(g_{i}^{\mathbf{t}})$
$\displaystyle=$ $\displaystyle\tau_{i}.$ (9)

The first equality follows from the definition of $Y_{i}$, the second equality
from the result in (3) and A3, and the third equality follows from
independence of $T$ and $Y$ given in A0:
$\mathbb{E}[\epsilon_{i}|T_{i}]=\mathbb{E}[\epsilon_{i}]=\alpha$ for all $i$.
Finally, we can use both of the results above to obtain the ADE:

$\displaystyle\frac{1}{n^{(1)}}\sum_{i=1}^{n}\mathbb{E}[T_{i}\times(\mathbb{E}[Y_{i}|G_{\mathcal{N}_{i}}^{\mathbf{T}}\simeq
g_{i}^{\mathbf{t}},T_{i}=1]$
$\displaystyle\qquad\qquad\quad-\mathbb{E}[Y_{i}|G_{\mathcal{N}_{i}}^{\mathbf{T}}\simeq
g_{i}^{\mathbf{t}},T_{i}=0])]$
$\displaystyle=\frac{1}{n^{(1)}}\sum_{i=1}^{n}\mathbb{E}[T_{i}\tau_{i}]$ (By
(9))
$\displaystyle=\frac{1}{n^{(1)}}\sum_{i=1}^{n}\mathbb{E}[T_{i}\times(\mathbb{E}[Y_{i}(1,\mathbf{0})-Y_{i}(0,\mathbf{0})])]$
(By (8))
$\displaystyle=\frac{1}{n^{(1)}}\sum_{i=1}^{n}\Pr(T_{i}=1)(\mathbb{E}[Y_{i}(1,\mathbf{0})-Y_{i}(0,\mathbf{0})])$
$\displaystyle=\frac{1}{n^{(1)}}\sum_{i=1}^{n}\frac{n^{(1)}}{n}(\mathbb{E}[Y_{i}(1,\mathbf{0})-Y_{i}(0,\mathbf{0})])$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[Y_{i}(1,\mathbf{0})-Y_{i}(0,\mathbf{0})],$

where $\Pr(T_{i}=1)=\frac{n^{(1)}}{n}$ by assumption of complete
randomization. ∎

###  B Additional Theoretical Results

We study the expected error for one AME match on subgraphs under two
assumptions: that the true weights for the AME objective (the weighted Hamming
distance) are known, and that the candidate units for matching all have
independently generated neighborhoods, and none of the units in these
neighborhoods are being matched. Additional information on this setting is
available in the proof.

######  Proposition 2.

(Oracle AME Error With Independent Neighborhoods) Suppose that there are $N$
independently generated graphs, each with $n$ vertices, and all i.i.d. from
the same distribution over graphs:
$\Pr(G_{1}=g_{1},\dots,G_{N}=g_{N})=\prod_{i=1}^{N}p(g_{i})$. Assume matches
are only allowed between units in different graphs. Suppose additionally that
$n^{(1)}$ randomly chosen units within each graph are assigned treatment, so
that $\Pr(\mathbf{T}_{i}=\mathbf{t}_{i})=\binom{n}{n^{(1)}}^{-1}$. Assume
further that interference functions obey the following: $|f(g)-f(h)|\leq
K\mathbf{w}^{T}\mathbb{I}[S(g)\neq S(h)]$, where $\mathbf{w}$ is a vector of
positive, real-valued, importance weights for the subgraphs counts, such that
$\|\mathbf{w}\|_{1}=M_{n}$ for some constant $0<M_{n}<\infty$, and such that
the condition above is satisfied for $\mathbf{w}$, and $\mathbb{I}[S(g)\neq
S(h)]$ is the Hamming distance: a vector of the same length as $\mathbf{w}$
that is 0 at position $k$ if graphs $g$ and $h$ have the same count of
subgraph $k$, and 1 otherwise. Assume that baseline responses have variance
$Var(\epsilon_{i})=\sigma^{2}\;\forall i$. Then, for a treatment unit $i$, if
$j$ solves the AME problem, i.e.,
$j\in\operatorname*{arg\,min}\limits_{\begin{subarray}{c}{k=1,\dots,n}\\\
{T_{k}=0}\end{subarray}}\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}})\neq
S(G_{\mathcal{N}_{k}}^{\mathbf{T}})]$, under A0-A3:

$\displaystyle\mathbb{E}\left[|Y_{i}-Y_{j}-\tau_{i}|\big{|}T_{i}=1,G_{\mathcal{N}_{i}}^{\mathbf{t}}=h_{\mathcal{N}_{i}}^{\mathbf{t}}\right]\leq\sqrt{2}\sigma$
$\displaystyle+K\binom{n-1}{n^{(1)}}^{-1}\times\sum_{g\in\mathcal{G}_{n}}\sum_{\begin{subarray}{c}\mathbf{t}\in\mathcal{T}\\\
t_{j}=0\end{subarray}}\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}})\neq
S(g_{\mathcal{N}_{j}}^{\mathbf{t}})]p(g)$
$\displaystyle\times\left(\frac{n^{(1)}}{n}+\frac{n-n^{(1)}}{n}C(g_{\mathcal{N}_{j}}^{\mathbf{t}})\right)^{N-2},$

where $\mathcal{G}_{n}$ is the set of all graphs with $n$ units, and
$C(h_{\mathcal{N}_{j}}^{\mathbf{t}})\leq 1$ for all $g$ and $\mathbf{t}$.

A proof is available in the following section. The first element in the right
hand side of the inequality is the standard deviation of the baseline
responses. One summation is over all possible graphs with $n$ units, and the
other summation is over possible treatment assignments. The expression inside
the summation is the product of three terms. First, the weighted Hamming
distance between a graph and the target graph we are trying to match. Second,
the probability of observing that graph. Third, an upper bound on the
probability that unit $j$ is among the minimizers of the weighted Hamming
distance. Note that
$\mathbf{w}^{T}\mathbb{I}[S(g_{\mathcal{N}_{j}}^{\mathbf{t}})\neq
S(g_{\mathcal{N}_{j}}^{\mathbf{t}})]p(g)$ is bounded for fixed $n$ for all $g$
and $\mathbf{t}$. This implies that the bound converges to $2\sqrt{\sigma}$ as
$N\rightarrow\infty$, as long as the size of neighborhood graphs is held
fixed, because perfect matching is possible with large amounts of data in this
regime.

####  B.1 Proof of Proposition 2

We briefly review some notation and assumptions to be used in this proof. For
the purposes of theory, we study a simplified setting, in which we have to
AME-match a unit $i$ to one unit in a set of candidate units of size $N$ such
that: a) all the candidate units belong to disconnected graphs, which we refer
to as candidate graphs. b) within each candidate graph there is only one pre-
determined candidate unit c) candidate units have neighborhood graphs denoted
by $G_{\mathcal{N}_{j}}$. d) all the candidate graphs are drawn independently
from the same distribution over graphs:
$\Pr(G_{\mathcal{N}_{1}}=g_{1},\dots,G_{\mathcal{N}_{N}}=g_{N})=\prod_{i=1}^{N}p(g_{i})$.
The support of $p$ will be $\mathcal{G}_{n}$: the set of all graphs with
exactly $n$ units. We use $g_{\mathcal{N}_{i}}$ to denote the subgraph induced
over $g$ by the units in the set of neighbors of unit $i$,
$\mathcal{N}_{i}\subseteq V(g)$, i.e., $g_{\mathcal{N}_{i}}$ is the graph
consisting only of the vertices that share an edge with $i$, in $g$, and of
the edges in $g$ that are between these vertices. The ego $i$ is not included
in $g_{\mathcal{N}_{i}}$.

Assigned treatments are denoted by $\mathbf{T}$, where
$\mathbf{T}\in\\{0,1\\}^{n}$, but in this setting treatment assignment is
assumed to be independent within the $N$ candidate graphs. Formally, the
assumption we make is that
$\Pr(\mathbf{T}_{1}=\mathbf{t}_{1},\dots,\mathbf{T}_{N}=\mathbf{t}_{N})=\prod_{i=1}^{N}\binom{n}{n^{(1)}}^{-1}$,
i.e., $n^{(1)}$ units are always treated uniformly at random within each of
the $N$ candidate graphs.

The direct treatment effect for any unit $i$ is given by $\tau_{i}$. We use
$\mathbb{I}[S(g)\neq S(h)]$ to indicate the Hamming distance between subgraph
counts of graphs $g$ and $h$. This means that $\mathbb{I}[S(g)\neq S(h)]$ is a
vector of size $|\mathcal{G}_{n}|$ that will be 1 in the $\ell^{th}$ entry if
$g$ and $h$ have the same amount of occurrences of graph $g_{\ell}$ among
their subgraphs. Note that this distance is coloring sensitive: two subgraphs
that are isomorphic in shape but not labels will belong to different entries
in this distance. The matched group of a treated unit $i$, denoted
$\mathtt{MG}_{i}$ is the set of all control units that match with $i$. In our
setting $j\in\mathtt{MG}_{i}$ if it solves the AME problem, that is
$j\in\operatorname*{arg\,min}\limits_{\begin{subarray}{c}k=1\dots,n\\\
T_{k}=0\end{subarray}}\mathbf{w}^{T}\mathbb{I}[S(G_{\mathcal{N}_{k}}^{\mathbf{T}})\neq
S(g_{\mathcal{N}_{i}}^{\mathbf{t}})]$. Finally, we assume that both the graph
for the unit we want to match and the treatment assignment for that unit’s
graph are fixed: $\mathbf{t}_{i}$ is the treatment assignment in the graph of
$i$, and $h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}}$ is the neighborhood graph of
$i$, where $h$ denotes unit $i$’s graph. All other notation is as in the main
paper.

###### Proof.

We start by upper-bounding our quantity of interest as follows:

$\displaystyle\mathbb{E}[|Y_{i}-Y_{j}-\tau_{i}|\big{|}j\in\mathtt{MG}_{i}]$
$\displaystyle=\mathbb{E}[|Y_{i}(1,\mathbf{t}_{i-i})-Y_{j}(0,\mathbf{T}_{j-j})-\tau_{i}|\big{|}j\in\mathtt{MG}_{i}]$
$\displaystyle=\mathbb{E}[|\tau_{i}+f(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})+\epsilon_{i}-f(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}})-\epsilon_{j}-\tau_{i}|\big{|}j\in\mathtt{MG}_{i}]$
$\displaystyle\leq\mathbb{E}[|f(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})-f(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}})|\big{|}j\in\mathtt{MG}_{i}]$
$\displaystyle\quad+\mathbb{E}[|\epsilon_{i}-\epsilon_{j}|\big{|}j\in\mathtt{MG}_{i}]$
$\displaystyle\leq
K\mathbb{E}[\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}})]\big{|}j\in\mathtt{MG}_{i}]$
$\displaystyle\quad+\mathbb{E}[|\epsilon_{i}-\epsilon_{j}|\big{|}j\in\mathtt{MG}_{i}],$
(10)

where the notation $\mathbf{T}_{j-j}$ denotes the treatment indicator for
candidate graph $j$ for all units except $j$. The first equality follows from
A0 since the event $j\in\mathtt{MG}_{i}$ implies that $T_{j}=0$, as only
control units are allowed in the matched groups. The second equality follows
from Proposition 1. The first inequality is an application of the triangle
inequality. The last line follows from the condition on the interference
functions. Consider the second term. We can use the Cauchy-Schwarz inequality
to construct a simple upper bound on it:

$\displaystyle\mathbb{E}[|\epsilon_{i}-\epsilon_{j}|\big{|}j\in\mathtt{MG}_{i}]=\mathbb{E}[|\epsilon_{i}-\epsilon_{j}|]$
$\displaystyle\leq\sqrt{\mathbb{E}[(\epsilon_{i}-\epsilon_{j})^{2}]}$
$\displaystyle=\sqrt{\mathbb{E}[\epsilon_{i}^{2}]+\mathbb{E}[\epsilon_{j}^{2}]-2\mathbb{E}[\epsilon_{i}]\mathbb{E}[\epsilon_{j}]}=\sqrt{2}\sigma$

where the last equality follows for the fact that the $\epsilon_{i}$ have mean
0 and are independent, with $Var(\epsilon_{i})=\sigma^{2}$ for all $i$.

Consider now the term
$\mathbb{E}[\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}})]|j\in\mathtt{MG}_{i}]$. To upper-bound
this, we write it out as follows using the definition of expectation:

$\displaystyle\mathbb{E}[\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}})]|j\in\mathtt{MG}_{i}]$
$\displaystyle=\sum_{g\in\mathcal{G}_{n}}\sum_{\mathbf{t}\in\mathcal{T},t_{j}=0}\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(g_{\mathcal{N}_{j}}^{\mathbf{t}_{j}})]$
$\displaystyle\qquad\qquad\times\Pr(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}}=g_{\mathcal{N}_{j}}^{\mathbf{t}},\mathbf{T}_{j}=\mathbf{t}|j\in\mathtt{MG}_{i}).$

We want to find an upper bound on
$\Pr(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}}=g_{\mathcal{N}_{j}}^{\mathbf{t}},\mathbf{T}_{j}=\mathbf{t}|j\in\mathtt{MG}_{i})$.
We start by writing this quantity out as a product of two probabilities:

$\displaystyle\Pr(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}}=g_{\mathcal{N}_{j}}^{\mathbf{t}},\mathbf{T}_{j}=\mathbf{t}|j\in\mathtt{MG}_{i})$
$\displaystyle=\Pr(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}}=g_{\mathcal{N}_{j}}^{\mathbf{t}}|j\in\mathtt{MG}_{i},\mathbf{T}_{j}=\mathbf{t})$
$\displaystyle\times\Pr(\mathbf{T}_{j}=\mathbf{t}|j\in\mathtt{MG}_{i})$
$\displaystyle=\Pr(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}}=g_{\mathcal{N}_{j}}^{\mathbf{t}}|j\in\mathtt{MG}_{i},\mathbf{T}_{j}=\mathbf{t})\binom{n-1}{n^{(1)}}^{-1}.$

Note that
$\Pr(\mathbf{T}_{j}=\mathbf{t}|j\in\mathtt{MG}_{i})=\binom{n-1}{n^{(1)}}^{-1}$
because treatment is uniformly randomized with $n^{(1)}$ units always treated
in each candidate graph, but $T_{j}=0$ conditionally on $j\in\mathtt{MG}_{i}$.

We use Bayes’ rule to write out the first term in the final product as
$\Pr(G_{\mathcal{N}_{j}}^{\mathbf{T}}=g_{\mathcal{N}_{j}}^{\mathbf{t}}|j\in\mathtt{MG}_{i},\mathbf{T}_{j}=\mathbf{t})=\frac{\Pr(j\in\mathtt{MG}_{i}|\mathbf{T}_{j}=\mathbf{t},G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}}=g_{\mathcal{N}_{j}}^{\mathbf{t}})\Pr(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}}=g_{\mathcal{N}_{j}}^{\mathbf{t}}|\mathbf{T}_{j}=\mathbf{t})}{\Pr(j\in\mathtt{MG}_{i}|\mathbf{T}_{j}=\mathbf{t})}$.

By assumption, if all neighborhood graphs are empty, all units are used for
all matched groups, and we are restricting ourselves to assignments in which
$T_{j}=0$, therefore, $\Pr(j\in\mathtt{MG}_{i}|\mathbf{T}_{j}=\mathbf{t})=1$.
Second, by assumption
$\Pr(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}}=g_{\mathcal{N}_{j}}^{\mathbf{t}}|\mathbf{T}_{j}=\mathbf{t})=p(g)$.
This is because treatment assignment is independent of the graph. We are left
with having to find an expression for the likelihood, this can be written as:

$\displaystyle\Pr(j\in\mathtt{MG}_{i}|\mathbf{T}_{j}=\mathbf{t},G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}}=g_{\mathcal{N}_{j}}^{\mathbf{t}})$
$\displaystyle=\Pr(j\in\operatorname*{arg\,min}_{\begin{subarray}{c}k=1,\dots,N,\\\
k\neq
i\end{subarray}}\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(G_{\mathcal{N}_{k}}^{\mathbf{T}_{k}})]$
$\displaystyle\qquad\qquad|\mathbf{T}_{j}=\mathbf{t},G_{\mathcal{N}_{j}}^{\mathbf{T}}=g_{\mathcal{N}_{j}}^{\mathbf{t}})$
$\displaystyle=\prod_{\begin{subarray}{c}k=1\\\ k\neq
i,j\end{subarray}}^{N}\bigg{[}\Pr(\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(G_{\mathcal{N}_{k}}^{\mathbf{T}_{k}})]\geq$
$\displaystyle\qquad\qquad\>\quad\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(g_{\mathcal{N}_{j}}^{\mathbf{t}})]|T_{k}=0)\Pr(T_{k}=0)$
$\displaystyle\qquad\quad+\Pr(T_{k}=1)\bigg{]}$
$\displaystyle=\prod_{\begin{subarray}{c}k=1\\\ k\neq
i,j\end{subarray}}^{N}\bigg{[}\Pr\big{(}\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(G_{\mathcal{N}_{k}}^{\mathbf{T}_{k}})]$ $\displaystyle\qquad\quad\ \
\,\geq\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(g_{\mathcal{N}_{j}}^{\mathbf{t}})]|T_{k}=0\big{)}\frac{n-n^{(1)}}{n}$
$\displaystyle\qquad\quad+\frac{n^{(1)}}{n}\biggl{]}.$

The second equality follows because $k$ can never be in the matched group of
unit $i$ if $T_{k}=1$, and, if $T_{k}=0$, then $k$ must be one of the
minimizers of the weighted Hamming distance between neighborhood subgraph
counts. The probability is a product of densities because of independence of
candidate subgraphs. For an arbitrary unit, $k$, we define the following
compact notation for the probability that $k$’s weighted Hamming distance from
$i$ is larger than the weighted Hamming distance from $j$ to $i$:

$\displaystyle\Pr(\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(G_{\mathcal{N}_{k}}^{\mathbf{T}_{k}})]$
$\displaystyle\qquad\geq\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(g_{\mathcal{N}_{j}}^{\mathbf{t}})]|T_{k}=0)$
$\displaystyle=:C_{k}(g_{\mathcal{N}_{j}}^{\mathbf{t}})\leq 1.$

Note that the last inequality follows from the fact that the expression above
is a probability. Since graphs and treatment assignments, $G_{k}$ and
$\mathbf{T}_{k}$ are the only random variables in the probability denoted by
$C_{k}(g_{\mathcal{N}_{j}}^{\mathbf{t}})$, and since they are all independent,
and identically distributed, we can say that
$C_{1}(g_{\mathcal{N}_{j}}^{\mathbf{t}})=C_{2}(g_{\mathcal{N}_{j}}^{\mathbf{t}})=\dots=C_{N}(g_{\mathcal{N}_{j}}^{\mathbf{t}})=C(g_{\mathcal{N}_{j}}^{\mathbf{t}})$.
Because of this we have:

$\displaystyle\Pr(j\in\mathtt{MG}_{i}|G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}}=g_{\mathcal{N}_{j}}^{\mathbf{t}},\mathbf{T}_{j}=\mathbf{t})$
$\displaystyle=\left(\frac{n^{(1)}}{n}+\frac{n-n^{(1)}}{n}C(g_{\mathcal{N}_{j}}^{\mathbf{t}})\right)^{N-2}.$

Putting all the elements we have together we get the expression for the first
term in the bound:

$\displaystyle\mathbb{E}[\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}})]|j\in\mathtt{MG}_{i}]$
$\displaystyle=\sum_{g\in\mathcal{G}_{n}}\sum_{\mathbf{t}\in\mathcal{T}:t_{j}=0}\mathbf{w}^{T}\mathbb{I}[S(g_{\mathcal{N}_{j}}^{\mathbf{t}})\neq
S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})]$
$\displaystyle\qquad\qquad\times\Pr(G_{\mathcal{N}_{j}}^{\mathbf{T}}=g_{\mathcal{N}_{j}}^{\mathbf{t}},\mathbf{T}_{\mathcal{N}_{j}}=\mathbf{t}_{\mathcal{N}_{j}}|j\in\mathtt{MG}_{i})$
$\displaystyle=\sum_{g\in\mathcal{G}_{N}}\sum_{\mathbf{t}\in\mathcal{T}:t_{j}=0}\mathbf{w}^{T}\mathbb{I}[S(g_{\mathcal{N}_{j}}^{\mathbf{t}})\neq
S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})]$
$\displaystyle\times\binom{n-1}{n^{(1)}}^{-1}p(g)\left(\frac{n^{(1)}}{n}+\frac{n-n^{(1)}}{n}C(g_{\mathcal{N}_{j}}^{\mathbf{t}})\right)^{N-2}.$

∎

####  B.2 Asymptotic Behavior

Here we expand on the asymptotic consequences of Proposition 2: note first,
that, by assumption $\|\mathbf{w}\|_{1}=M_{n}$, and that, therefore
$\mathbf{w}^{T}\mathbb{I}[S(g)\neq S(h)]\leq M_{n}$ for any graphs
$g,h\in\mathcal{G}_{n}$. That is to say, the weighted Hamming distance between
any two graphs with $n$ units will be upper-bounded by the sum of the weights.
Recall also that $C(g_{\mathcal{N}_{j}}^{\mathbf{t}})\leq 1$ for all $g$ and
$\mathbf{t}$ as this quantity is a probability, and let
$C_{max}=\max\limits_{\begin{subarray}{c}g\in\mathcal{G}_{n}\\\
\mathbf{t}\in\mathcal{T}\end{subarray}}C(g_{\mathcal{N}_{j}}^{\mathbf{t}})$.
We can combine all these bounds with the upper bound in Proposition 2 to
write:

$\displaystyle\mathbb{E}[\mathbf{w}^{T}\mathbb{I}[S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})\neq
S(G_{\mathcal{N}_{j}}^{\mathbf{T}_{j}})]|j\in\mathtt{MG}_{i}]$
$\displaystyle=\sum_{g\in\mathcal{G}_{\mathcal{N}}}\sum_{\mathbf{t}\in\mathcal{T},t_{j}=0}\mathbf{w}^{T}\mathbb{I}[S(g_{\mathcal{N}_{j}}^{\mathbf{t}})\neq
S(h_{\mathcal{N}_{i}}^{\mathbf{t}_{i}})]$
$\displaystyle\times\binom{n-1}{n^{(1)}}^{-1}p(g)\left(\frac{n^{(1)}}{n}+\frac{n-n^{(1)}}{n}C(g_{\mathcal{N}_{j}}^{\mathbf{t}})\right)^{N-2}$
$\displaystyle\leq
M_{n}\left(\frac{n^{(1)}}{n}+\frac{n-n^{(1)}}{n}C_{max}\right)^{N-2}$
$\displaystyle\times\sum_{g\in\mathcal{G}_{n}}\sum_{\mathbf{t}\in\mathcal{T},t_{j}=0}\binom{n-1}{n^{(1)}}^{-1}p(g)$
$\displaystyle=M_{n}\left(\frac{n^{(1)}}{n}+\frac{n-n^{(1)}}{n}C_{max}\right)^{N-2}.$

The first equality follows from Proposition 2, the first inequality from the
bounds previously discussed, and the second equality follows from the fact
that the sum in the second to last line is a sum of probability distributions
over their entire domain, and therefore is equal to 1. Under the condition
that $n$, the number of units in each unit’s candidate graph, stays fixed, and
that $C_{max}<1$, then, as $N\rightarrow\infty$, we have
$M_{n}\left(\frac{n^{(1)}}{n}+\frac{n-n^{(1)}}{n}C_{max}\right)^{N-2}\rightarrow
0$, because the quantity inside the parentheses is always less than 1. This
makes sense, because asymptotically, matches can be made exactly; i.e., units
matched in the way described in our theoretical setting have isomorphic
neighborhood subgraphs asymptotically. This also has a consequence that the
bound in Proposition 2 converges to $\sqrt{2}\sigma$ asymptotically in $N$.
This is the variance of the baseline errors and can be lowered by matching the
same unit with multiple others. As noted before, for this argument to apply,
candidate graphs must remain of fixed size $n$ as they grow in number, so that
the quantity $M_{n}$ remains constant: this setting is common in cluster-
randomized trials where a growing number of units is clustered into fixed-size
clusters of size at most $n$. The asymptotic behavior of our proposed
methodology is less clear in settings in which the analyst cannot perform such
clustering before randomization and $n$ is allowed to grow with $N$, and is an
avenue for potential future research.

####  B.3 Heteroskedasticity in The Baseline Effects

In a network setting such as the one we study, it is possible that baseline
effects of units do not have equal variance. Here we discuss how this setting
affects our result in Proposition 2. Here, we maintain that
$\mathbb{E}[\epsilon_{i}]=\alpha$ for all $i$, but we assume that
$Var(\epsilon_{i})=\sigma_{i}^{2}$, and that
$Cov(\epsilon_{i}\epsilon_{j})\neq 0$. Starting from the upper bound on the
estimation error given in (10), we can see that the baseline effects only come
in in the term: $\mathbb{E}[|\epsilon_{i}-\epsilon_{j}||j\in\mathtt{MG}_{i}]$,
we therefore focus our attention on this term, as the rest of this bound does
not change when the variance of these terms changes. Note first, that
$\mathbb{E}[|\epsilon_{i}-\epsilon_{j}||j\in\mathtt{MG}_{i}]=\mathbb{E}[|\epsilon_{i}-\epsilon_{j}|]$
as the event $j\in\mathtt{MG}_{i}$ is independent of the baseline effects. We
can now apply the Cauchy-Schwarz inequality, in the same way as we do in the
proof of Proposition 2, to obtain:

$\displaystyle\mathbb{E}[|\epsilon_{i}-\epsilon_{j}|]$
$\displaystyle\leq\sqrt{\mathbb{E}[(\epsilon_{i}-\epsilon_{j})^{2}]}$
$\displaystyle=\sqrt{\mathbb{E}[\epsilon_{i}^{2}]+\mathbb{E}[\epsilon_{j}^{2}]-2\mathbb{E}[\epsilon_{i}]\mathbb{E}[\epsilon_{j}]}$
$\displaystyle=\sqrt{\sigma^{2}_{i}+\alpha^{2}+\sigma_{j}^{2}+\alpha^{2}-2\alpha^{2}}$
$\displaystyle=\sqrt{\sigma_{i}^{2}+\sigma_{j}^{2}}.$

Clearly, this is not too different from the homoskedastic setting we study in
the proposition: as long as neither of the unit variances is too large for
inference, results in the heteroskedastic setting will suffer from similar
bias as they would under independent baseline effects with equal variance.

Simulations, shown in Section J, also support the above rationale of
comparable performance in the heteroskedastic case and demonstrate that FLAME-
Networks still outperforms competing methods.

###  C Derivation of SANIA MIVLUE Used in Simulations

Theorem 6.2 of Sussman and Airoldi, 2017
Suppose potential outcomes satisfy SANIA and that the prior on the parameters
(baseline outcome and direct treatment effect) has no correlation between
units. If unbiased estimators exist, the MIVLUE weights are:

$w_{i}(\mathbf{z})=\frac{C_{i,d_{i}^{\mathbf{z}}}}{\sum_{d=0}^{n-1}C_{i,d}}\cdot\frac{2z_{i}-1}{n\mathrm{P}(\mathbf{z}_{i}^{\text{obs}}=z_{i},d_{i}^{\mathbf{z}^{\text{obs}}}=d_{i}^{\mathbf{z}})}$

where

$\displaystyle C_{i,d}$ $\displaystyle=$
$\displaystyle\left(\sum_{\mathbf{z}}\mathrm{P}(\mathbf{z})\mathbf{1}\\{d_{i}^{\mathbf{z}}=d\\}\cdot\frac{\Sigma(\mathbf{z})_{ii}}{\mathrm{P}(z_{i}^{\text{obs}}=z_{i},d_{i}^{\mathbf{z}^{\text{obs}}}=d)^{2}}\right)^{-1}$

and $C_{i,d}$ is defined to be 0 if the probability in its denominator is 0.

####  C.1 Setup and Notation

We assume that there is a constant probability $p$ of each unit being treated
and that units are independently treated. Let unit $i$ have $d_{i}$ neighbors
and a treated degree of $d_{i}^{\mathbf{z}}$. In our setting,
$\Sigma(\mathbf{z})_{ii}=\Sigma_{\alpha,ii}+z_{i}\Sigma_{\beta,ii}$ where
$\Sigma_{\alpha}$ and $\Sigma_{\beta}$ are the covariance matrices on priors
placed on the baseline outcome and the direct treatment effect, respectively.
Additionally in our setting, their diagonals are constant and so we let
$\sigma^{2}_{\alpha}:=\Sigma_{\alpha,ii}$ and
$\sigma^{2}_{\beta}:=\Sigma_{\beta,ii}$.

####  C.2 Find
$\mathrm{P}(z_{i}^{\text{obs}}=z_{i},d_{i}^{\mathbf{z}^{\text{obs}}}=d_{i}^{\mathbf{z}})$

By the setup, the constituent probabilities are independent and the
probability of treatment is constant across units and so:
$\mathrm{P}(z_{i}^{\text{obs}}=z_{i},d_{i}^{\mathbf{z}^{\text{obs}}}=d_{i}^{\mathbf{z}})=[z_{i}p+(1-z_{i})(1-p)]{d_{i}\choose
d_{i}^{\mathbf{z}}}p^{d_{i}^{\mathbf{z}}}(1-p)^{d_{i}-d_{i}^{\mathbf{z}}}$

####  C.3 Find $C_{i,d}$

Below, we will consider $\mathbf{z}_{\text{neighbor}(i)}$, the treatment
assignment of the neighbors of $i$ (excluding $i$), $z_{i}$, the treatment
assignment of unit $i$, and $\mathbf{z}_{\text{rest}(i)}$, the treatment
assignment of the remaining units.

$\displaystyle C_{i,d}$
$\displaystyle=\left(\sum_{\mathbf{z}}\frac{\mathrm{P}(\mathbf{z})\mathbf{1}\\{d_{i}^{\mathbf{z}}=d\\}\Sigma(\mathbf{z})_{ii}}{\mathrm{P}(z_{i}^{\text{obs}}=z_{i},d_{i}^{\mathbf{z}^{\text{obs}}}=d)^{2}}\right)^{-1}$
$\displaystyle=\left(\sum_{\mathbf{z}:d_{i}^{\mathbf{z}}=d}\frac{\mathrm{P}(\mathbf{z}_{\text{neighbor}(i)})\mathrm{P}(z_{i})\mathrm{P}(\mathbf{z}_{\text{rest}(i)})\Sigma(\mathbf{z})_{ii}}{\mathrm{P}(z_{i}^{\text{obs}}=z_{i},d_{i}^{\mathbf{z}^{\text{obs}}}=d)^{2}}\right)^{-1}$
$\displaystyle=\left(\frac{\mathrm{P}(\mathbf{z}_{\text{neighbor}(i)})}{\mathrm{P}(z_{i}^{\text{obs}}=z_{i},d_{i}^{\mathbf{z}^{\text{obs}}}=d)^{2}}\right)^{-1}$
$\displaystyle\quad\times\left(\sum_{\mathbf{z}:d_{i}^{\mathbf{z}}=d}\Sigma(\mathbf{z})_{ii}\mathrm{P}(z_{i})\mathrm{P}(\mathbf{z}_{\text{rest}(i)})\right)^{-1}$
$\displaystyle=\left(\frac{\mathrm{P}(\mathbf{z}_{\text{neighbor}(i)})}{\mathrm{P}(z_{i}^{\text{obs}}=z_{i},d_{i}^{\mathbf{z}^{\text{obs}}}=d)^{2}}\right)^{-1}$
$\displaystyle\quad\times\left(\sum_{\mathbf{z}:d_{i}^{\mathbf{z}}=d}(\sigma^{2}_{\alpha}+z_{i}\sigma^{2}_{\beta})\mathrm{P}(z_{i})\mathrm{P}(\mathbf{z}_{\text{rest}(i)})\right)^{-1}$
$\displaystyle=\left(\frac{\mathrm{P}(\mathbf{z}_{\text{neighbor}(i)})}{\mathrm{P}(z_{i}^{\text{obs}}=z_{i},d_{i}^{\mathbf{z}^{\text{obs}}}=d)^{2}}\right)^{-1}$
$\displaystyle\quad\times\left(\sigma^{2}_{\alpha}+\sum_{\mathbf{z}:d_{i}^{\mathbf{z}}=d}(z_{i}\sigma^{2}_{\beta})\mathrm{P}(z_{i})\mathrm{P}(\mathbf{z}_{\text{rest}(i)})\right)^{-1}$
$\displaystyle=\left(\frac{\mathrm{P}(\mathbf{z}_{\text{neighbor}(i)})}{\mathrm{P}(z_{i}^{\text{obs}}=z_{i},d_{i}^{\mathbf{z}^{\text{obs}}}=d)^{2}}\right)$
$\displaystyle\quad\times\left(\sigma^{2}{\alpha}+\sum_{\begin{subarray}{c}\mathbf{z}:d_{i}^{\mathbf{z}}=d\\\
z_{i}=1\end{subarray}}(\sigma^{2}_{\beta,})p\mathrm{P}(\mathbf{z}_{\text{rest}(i)})\right)^{-1}$
$\displaystyle=\left(\frac{\mathrm{P}(\mathbf{z}_{\text{neighbor}(i)})\left(\sigma^{2}_{\alpha}+\sigma^{2}_{\beta}\right)}{\mathrm{P}(z_{i}^{\text{obs}}=z_{i},d_{i}^{\mathbf{z}^{\text{obs}}}=d)^{2}}\right)^{-1}$
$\displaystyle=\frac{\left([z_{i}p+(1-z_{i})(1-p)]{d_{i}\choose
d}p^{d}(1-p)^{d_{i}-d}\right)^{2}}{(\sigma^{2}_{\alpha}+\sigma^{2}_{\beta})p^{d}(1-p)^{d_{i}-d}}$
$\displaystyle=\frac{[z_{i}p+(1-z_{i})(1-p)]^{2}{d_{i}\choose
d}^{2}p^{d}(1-p)^{d_{i}-d}}{\sigma^{2}_{\alpha}+\sigma^{2}_{\beta}}$

####  C.4 Find $w_{i}$

Plugging in the expressions we’ve found:

$\displaystyle w_{i}(\mathbf{z})=$
$\displaystyle\frac{C_{i,d_{i}^{\mathbf{z}}}}{\sum_{d=0}^{n-1}C_{i,d}}\cdot\frac{2z_{i}-1}{n\mathrm{P}(\mathbf{z}_{i}^{\text{obs}}=z_{i},d_{i}^{\mathbf{z}^{\text{obs}}}=d_{i}^{\mathbf{z}})}$
$\displaystyle=\frac{\frac{[z_{i}p+(1-z_{i})(1-p)]^{2}{d_{i}\choose
d_{i}^{\mathbf{z}}}^{2}p^{d_{i}^{\mathbf{z}}}(1-p)^{d_{i}-d_{i}^{\mathbf{z}}}}{(\sigma^{2}_{\alpha}+\sigma^{2}_{\beta})}}{\frac{\sum_{d=0}^{n-1}{d_{i}\choose
d}^{2}p^{d}(1-p)^{d_{i}-d}(z_{i}p+(1-z_{i})(1-p))^{2}}{(\sigma^{2}_{\alpha}+\sigma^{2}_{\beta})}}$
$\displaystyle\quad\times\frac{2z_{i}-1}{n[z_{i}p+(1-z_{i})(1-p)]{d_{i}\choose
d_{i}^{\mathbf{z}}}p^{d_{i}^{\mathbf{z}}}(1-p)^{d_{i}-d_{i}^{\mathbf{z}}}}$
$\displaystyle=\frac{(z_{i}p+(1-z_{i})(1-p))(2z_{i}-1)}{n\sum_{d=0}^{n-1}{d_{i}\choose
d}^{2}p^{d}(1-p)^{d_{i}-d}(z_{i}p+(1-z_{i})(1-p))^{2}}$
$\displaystyle=\frac{(z_{i}p+(1-z_{i})(1-p))(2z_{i}-1)}{n(z_{i}p+(1-z_{i})(1-p))^{2}}$
$\displaystyle\quad\times\frac{1}{\sum_{d=0}^{n-1}{d_{i}\choose
d}^{2}p^{d}(1-p)^{d_{i}-d}}$

Note in the first fraction that $(z_{i}p+(1-z_{i})(1-p))(2z_{i}-1)$ equals $p$
when $z_{i}=1$ and $p-1$ when $z_{i}=0$. Also, $(z_{i}p+(1-z_{i})(1-p))^{2}$
equals $p^{2}$ when $z_{i}=1$ and $(1-p)^{2}$ when $z_{i}=0$. Thus, the first
term is $1/np$ when $z_{i}=1$ and $-1/n(1-p)$ when $z_{i}=0$ and so the
overall expression for the weights is:

$w_{i}(\mathbf{z})=\frac{z_{i}/np-(1-z_{i})/(n(1-p))}{\sum_{d=0}^{n-1}{d_{i}\choose
d}^{2}p^{d}(1-p)^{d_{i}-d}}$

and the MIVLUE is given by $\sum_{i=1}^{n}w_{i}Y_{i}$.

###  D Subgraph Descriptions

Here, for graphs without self-loops, we define the interference components
used in our simulations:

  * •

Degree: the degree of a node is the number of edges it is a part of.

  * •

Triangle: A graph with three mutually connected nodes (see Figure 5).

  * •

Square: A graph with 4 nodes and 4 edges, such that each node is a part of
exactly two distinct edges (see Figure 5).

  * •

$k$-Star: A graph with $k+1$ nodes, the first $k$ of which are all connected
to the $(k+1)$st node and no others (see Figure 5).

  * •

Vertex Betweenness: The vertex betweenness of a vertex $v$ is defined as:

$\sum_{i\neq v\neq j}\frac{\sigma_{ij}(v)}{\sigma_{ij}}$

where $\sigma_{ij}$ is the number of shortest paths between vertices $i$ and
$j$ and $\sigma_{ij}(v)$ is the number of shortest paths between $i$ and $j$
that go through $v$. We use the normalized vertex betweenness which scales the
above expression by $2/(n^{2}-3n+2)$ where $n$ is the number of nodes in the
graph.

  * •

Closeness Centrality: We use the normalized closeness centrality of a vertex
$v$, defined as:

$\frac{n-1}{\sum_{i=1}^{n}d(\sigma_{vi})}$

where $d(\sigma_{vi})$ is the length of the shortest path between $v$ and $i$
and $n$ is the number of nodes in the graph.

Figure 5: Four types of treated neighborhood subgraphs the colored unit might
be a part of: triangle (top right), square (bottom right), 2-star (top left),
4-star (bottom left).

###  E Data Pre-processing

We estimate the DTE using data from from Banerjee et al. (2013) on social
networks for the 75 villages in India. A unit $i$ is defined to be socially
connected with unit $j$ if they are connected across at least 3 of the
following four types of social connections: (1) $i$ visits $j$’s house, (2)
$j$ visits $i$’s house, (3) $i$ and $j$ socialize and are relatives, (4) $i$
and $j$ socialize and are not related to each other. Along with subgraph
counts, we also use age and whether or not individual $i$ speaks Telugu as
additional covariates to match on. For computational efficiency, we drop units
with maximum degree of connection greater than 15, where the cut-off is
selected based on computational efficiency.

###  F Matched Groups

We provide sample matched groups in Table 4. These matched groups were
produced by applying FLAME-Networks on the data discussed in Section 4. We
report all the covariates used for matching. The first group is comprised of
40-year-old units who do not speak Telugu, and have 2 or 3 triangles, 3
2-stars, and 3 edges in their treated neighborhood graph. These units (given
the binning) are matched exactly. The second group is comprised of units who
speak Telugu with no triangles, 3 2-stars, 3 edges in their treated
neighborhood graph. Note that units in this group are matched approximately,
since they are not matched exactly on age.

Units | Triangles | 2-Stars | Edges | Telugu | Age | Treated | Outcome
---|---|---|---|---|---|---|---
Matched Group 1 |  |  |  |  |  |  |
1 | 2 or 3 | 3 | 3 | 0 | 40 | 0 | 1
2 | 2 or 3 | 3 | 3 | 0 | 40 | 1 | 0
3 | 2 or 3 | 3 | 3 | 0 | 40 | 1 | 0
4 | 2 or 3 | 3 | 3 | 0 | 40 | 1 | 0
Matched Group 2 |  |  |  |  |  |  |
1 | 0 | 3 | 3 | 1 | 30 | 0 | 0
2 | 0 | 3 | 3 | 1 | 34 | 0 | 0
3 | 0 | 3 | 3 | 1 | 25 | 0 | 1
4 | 0 | 3 | 3 | 1 | 20 | 1 | 0

Table 4: Sample Match Groups. Two sample matched groups generated by FLAME-
Networks using data discussed in Section 4. The columns are the covariates
used for matching, along with treatment status and outcome. The counts of
subgraphs were coarsened into 10 bins defined by deciles of the empirical
distribution of counts. The two groups have relatively good match quality
overall. Note that the first group matches units exactly (given the binning).
However, Matched Group 2 matches units approximately, with exact matches on
subgraph counts and whether or not individuals speak Telugu, but inexact
matches on age.

###  G Additional Experiment: Multiplicative Interference

We explore settings in which interference is a nonlinear function of the
interference components and their weights. Since matching is nonparametric, it
is particularly appropriate for handling non-linearities in interference
functions. Outcomes in this experiment have the form
$Y_{i}(t,\mathbf{t}_{-i})=t\tau_{i}+\alpha\prod_{p=1}^{P}m_{ip}^{\mathbb{I}[p\text{
is included}]}+\epsilon_{i}$. Table 5 shows which components are included in
the outcome function for each setting. We use a small number of parameters in
each setting, as their multiplicative interaction suffices to define a complex
interference function. The simulations are run on an $ER(50,0.05)$ graph.

Setting | $d_{i}$ | $\Delta_{i}$ | $\bigstar_{i}^{4}$ | $B_{i}$
---|---|---|---|---
1 | x | x |  |
2 | x |  |  | x
3 |  | x |  | x
4 |  | x | x |

Table 5: Parameters included in interference function Experiment 1. The marked
components for each setting were the only ones included in those experiments.

Results for this experiment are presented in Figure 6. FLAME-Networks
performed better than all baseline methods in this setting, both in terms of
mean absolute error and, in most cases, in terms of standard deviation over
simulations. The stratified and SANIA estimators perform especially poorly,
because they cannot handle nonlinear interference settings, unlike FLAME-
Networks.

Figure 6: Results from experiments with a multiplicative interference
function. Each violin plot represents the distribution over simulations of
absolute estimation error over for each method. The panels are numbered
according to the parameter setting the simulations were ran with. Violin plots
are color-coded blue if the method had mean error either equal to or better
than FLAME-Networks and red otherwise. The black line inside each violin is
the median error. The dashed line is FLAME-Networks’ mean error.

###  H Additional Experiment: Graph Cluster Randomization

We also explored the performance of FLAME-Networks in settings in which
treatment is randomized within multiple clusters, which have few connections
between them. More specifically, we simulate a network according to a
stochastic block model with 5 clusters. In each cluster, there are 10 units, 5
of which are treated. The probability of edges within clusters is 0.3 and
between clusters is 0.05. This results in graphs with few edges between
clusters. We then evaluate our method as previously described, simulating the
outcome with additive interference and homoskedastic errors. The results in
Figure 7 demonstrate that FLAME-Networks outperforms competing methods in this
setting as well.

Figure 7: Results from experiments with additive interference on graphs in
which treatment is randomly assigned within multiple clusters with few edges
between them. Each violin plot represents the distribution over simulations of
absolute estimation error over for each method. The panels are numbered
according to the parameter setting the simulations were ran with. Violin plots
are color-coded blue if the method had mean error either equal to or better
than FLAME-Networks and red otherwise. The black line inside each violin is
the median error. The dashed line is FLAME-Networks’ mean error.

###  I Additional Experiment: Real Network

To ensure that FLAME-Networks also performs well on real networks, we consider
an AddHealth network (Harris, 2013). Specifically, we use the addhealthc3
dataset from the amen R package, with all edges treated as undirected. There
are 32 nodes and on every simulation, 16 are randomly selected to be treated.
Outcome and additive interference are simulated as previously described.
Errors are homoskedastic. The results in Figure 8 demonstrate that FLAME-
Networks still outperforms competing methods.

Figure 8: Results from experiments on a real, AddHealth network with additive
interference. Each violin plot represents the distribution over simulations of
absolute estimation error over for each method. The panels are numbered
according to the parameter setting the simulations were ran with. Violin plots
are color-coded blue if the method had mean error either equal to or better
than FLAME-Networks and red otherwise. The black line inside each violin is
the median error. The dashed line is FLAME-Networks’ mean error.

###  J Additional Experiment: Heteroscedastic Errors

We also explored the performance of FLAME-Networks in settings in which the
variance of the outcomes is not constant. We simulate a single ER(50, 0.07)
graph and randomly treat 25 units. We consider additive interference, as in
the body of the text, and all other simulation parameters are the same, expect
for that now, each unit $i$ has baseline outcome
$\alpha_{i}\stackrel{{\scriptstyle ind}}{{\sim}}N(0,v_{i})$ with
$v_{i}\stackrel{{\scriptstyle ind}}{{\sim}}U(0,1)$. We see in Figure 9 that
FLAME-Networks outperforms competitors; the fact that it is nonparametric
allows it to handle more flexible baseline outcomes and variances.

Figure 9: Results from experiments with an additive interference function
involving heteroskedasticity in the baseline effects across units. Each violin
plot represents the distribution over simulations of absolute estimation error
over for each method. The panels are numbered according to the parameter
setting the simulations were ran with. Violin plots are color-coded blue if
the method had mean error either equal to or better than FLAME-Networks and
red otherwise. The black line inside each violin is the median error. The
dashed line is FLAME-Networks’ mean error.

###  K Additional Experiment: Matching on True Interference

Here, we compare FLAME-Networks to approaches that match directly on units’
interference values. FLAME-Networks already has the advantage of interpretably
matching on neighborhood graphs that can be visualized and compared, as
opposed to uninterpretable scalar values of an interference function.
Additionally, to perform well in practice, one would typically need to use
equally uninterpretable machine learning methods to estimate units’
interference values well. But even comparing FLAME-Networks to an approach
that matches on the true (typically unknown) interference values, we see that
our method does well in comparison, because it learns and matches on baseline
effects as well as (approximate) interference values. Results using an ER(50,
0.07) graph with 25 units randomly treated, an additive interference function,
and homoskedastic errors – as previously described – are shown in Figure 10.

###  L Additional Experiment: Covariate Weight

In this section, we show that increasing the influence that covariates have on
the outcome function harms neither FLAME-Networks nor the competing methods.
As in the results shown in the main text, however, the performance of FLAME-
Networks is still superior, given that it naturally handles covariate data.
The experimental setup is the same as in Experiment 2 in the main text and
results are shown in Tables 6 and 7

Method | Median | 25th q | 75th q
---|---|---|---
FLAME-Networks | 0.34 | 0.15 | 0.52
First Eigenvector | 0.41 | 0.24 | 0.49
All Eigenvectors | 0.36 | 0.32 | 0.74
Naive | 0.61 | 0.19 | 0.85
SANIA | 2.31 | 1.78 | 2.75
Stratified | 4.56 | 4.55 | 4.63

Table 6: Additional results from the experimental setup of Experiment 2, but
with $\beta=20$. Median and 25th and 75th percentile of absolute error over 10
simulations.

–

Method | Median | 25th q | 75th q
---|---|---|---
FLAME-Networks | 0.25 | 0.08 | 0.51
First Eigenvector | 0.52 | 0.32 | 0.85
All Eigenvectors | 0.53 | 0.29 | 0.83
Naive | 0.78 | 0.32 | 1.16
SANIA | 1.86 | 1.68 | 2.11
Stratified | 4.78 | 4.74 | 4.80

Table 7: Additional results from the experimental setup of Experiment 2, but
with $\beta=25$. Median and 25th and 75th percentile of absolute error over 10
simulations.

–

Figure 10: Results from experiments comparing FLAME-Networks to matching units
on their true interference values. Each violin plot represents the
distribution over simulations of absolute estimation error over for each
method. The panels are numbered according to the parameter setting the
simulations were ran with. Violin plots are color-coded blue if the method had
mean error either equal to or better than FLAME-Networks and red otherwise.
The black line inside each violin is the median error. The dashed line is
FLAME-Networks’ mean error.

###  M Match Quality

Here, we assess the quality of matches generated by FLAME-Networks versus
matching on true interference, and the All Eigenvectors approach. To do so,
for FLAME-Networks: for each (control) treated unit, we take the minimal
Frobenius norm of the difference between that unit’s neighborhood adjacency
matrix and that of all the (treated) control units in its matched group111The
Frobenius norm of the difference of the adjacency matrices, up to reordering
of the vertices., and average across all units. This gives an average graph
distance for a single simulation. And to do so for the true interference
matching and All Eigenvectors approaches: for every (control) treated unit, we
take the closest (treated) control unit, find the graph distance (as above)
between their neighborhood subgraphs, and average across units. This gives an
average graph distance for a single simulation. Results from 50 simulations
performed on an ER(50, 0.07) graph with additive interference and
homoskedastic errors, as described previously, are shown in Figure 11, showing
that FLAME-Networks produces more matches between units with similar
neighborhood subgraphs than matching on the true interference or the All
Eigenvectors method.

Figure 11: Results from experiments comparing the average distance between the
neighborhood subgraphs of the units matched by different methods. Each violin
plot represents the distribution over simulations of graph distance for each
method. The panels are numbered according to the parameter setting the
simulations were ran with. Violin plots are color-coded blue if the method had
mean graph distance either equal to or better than FLAME-Networks and red
otherwise. The black line inside each violin is the median graph distance. The
dashed line is FLAME-Networks’ mean graph distance.