Causal inference for social network data

arXiv:1705.08527v4 [stat.ME] 5 Oct 2017

Elizabeth L. Ogburn∗, Oleg Sofrygin†, Iván Díaz‡, and Mark J. van der Laan†

Abstract
We extend recent work by van der Laan (2014) on causal inference for causally connected units
to more general social network settings. Our asymptotic results allow for dependence of each
observation on a growing number of other units as sample size increases. We are not aware of any
previous methods for inference about network members in observational settings that allow the
number of ties per node to increase as the network grows. While previous methods have generally
implicitly focused on one of two possible sources of dependence among social network observations,
we allow for both dependence due to contagion, or transmission of information across network ties,
and for dependence due to latent similarities among nodes sharing ties. We describe estimation
and inference for causal effects that are specifically of interest in social network settings.

1

Introduction and Background

Social networks have long been of interest to sociologists (Moreno, 1937) and in the past decade
they have gained prominence across a broad range of disciplines, from economics to epidemiology
(Knoke and Yang, 2008). A social network is, roughly, a set of individuals who share some kind
of social relationships with one another – e.g. friends, neighbors, family members, or coworkers. If
one individual’s outcome can be affected by his or her social contacts’ outcomes, then we say that the
outcome exhibits induction or contagion, and the causal effects are called peer effects. If an individual’s
outcome may be affected by his or her contacts’ treatments of exposures, then those treatments or
exposures are said to exhibit interference.
Many aspects of social networks are of interest to researchers, from the clustering of individuals
into communities to the probability distributions that describe the generation of new relationships
between individuals in the network. There is increasing interest in identifying and estimating causal
effects in the contexts of social networks, that is causal effects that one individual’s behavior, treatment
assignment, beliefs, or health outcome could have on his or her social contacts’ behaviors, exposures,
beliefs, or health statuses. There have been a number of high profile articles that use standard methods
like generalized linear models (GLM) and generalized estimating equations (GEE) to attempt to infer
causal peer effects from network data (e.g. Christakis and Fowler, 2007, 2008, 2010), and this work
∗ Department

of Biostatistics, Johns Hopkins Bloomberg School of Public Health, Baltimore, MD, USA
of Biostatistics, University of California, Berkeley, CA, USA
‡ Division of Biostatistics and Epidemiology, Weill Cornell Medicine, New York, NY, USA
† Department

1

has inspired several research programs that study peer effects using the same statistical methods (Ali
and Dwyer, 2010; Cacioppo et al., 2009; Madan et al., 2010; Rosenquist et al., 2010; Wasserman,
2013). However, these methods have come under considerable criticism from the statistical community
(Cohen-Cole and Fletcher, 2008; Lyons, 2011; Shalizi and Thomas, 2011). These statistical models
are not equipped to deal with dependence across individuals and are rarely appropriate for estimating
effects using network data (Ogburn and VanderWeele, 2014). In some settings it may be possible to use
them to test for the presence of network dependence, but some properties of such tests are unknown
(VanderWeele et al., 2012; Shalizi, 2012; Ogburn and VanderWeele, 2014).
Spatial autoregressive (SAR) models have been applied to the study of peer effects and induction in
network settings (e.g. Goetzke, 2008; Lee, 2004; Lin, 2005; O’Malley and Marsden, 2008). Because the
endogenous and exogenous variables are measured at the same time, they parameterize an equilibrium
state rather than causal relationships. Few data generating processes give rise to true equilibrium
states (Besag, 1974; Lauritzen and Richardson, 2002; Thomas, 2013); therefore SAR models may often
be misspecified or uninformative about causal relationships.
Very recently, researchers interested in causal inference for interconnected subjects have begun to
develop methods designed specifically for the network setting. Many methods for interference–the effect
of one individual’s treatment or exposure on others’ outcomes–are relevant to the analysis of network
data (Aronow and Samii, 2013; Athey et al., 2016; Bowers et al., 2013; Eckles et al., 2014; Graham
et al., 2010; Halloran and Struchiner, 1995; Halloran and Hudgens, 2011; Hong and Raudenbush, 2006,
2008; Hudgens and Halloran, 2008; Rosenbaum, 2007; Rubin, 1990; Sobel, 2006; Tchetgen Tchetgen and
VanderWeele, 2012; VanderWeele, 2010; VanderWeele and Tchetgen Tchetgen, 2011a,b). However, the
inferential methods developed in this context generally require observing multiple independent groups
of units, which corresponds to observing multiple independent networks, or else they require treatment
to be randomized. Ideally, we would like to be able to perform inference even when all observations
are sampled from a single social network and in observational settings in addition to randomized
experiments.
Methodology has not kept apace with interest in causal inference using data from individuals
connected in a social network, and researchers have resorted to using standard statistical methods to
analyze this new type of data, or to collecting multiple independent groups of individuals and using
the groups, rather than the individuals, as the unit of statistical inference. The former strategy is not
statistically appropriate and the the latter strategy, while statistically valid, can be an inefficient use
of data. Furthermore, it limits the settings and effects that one can study, because in some settings it
may be too expensive or labor intensive to collect data from many independent groups or independent
groups may not exist (e.g. a global infectious disease epidemic).
In this paper we extend recent work by van der Laan (2014) on causal inference for causally
connected units to more general social network settings. van der Laan (2014) introduced methods
for causal inference from a single collection of interconnected units when each unit is known to be
independent of all but a small number of other units. Asymptotic results rely on the number of
dependent units being fixed as the total number of units goes to infinity. We propose new methods
that allow similar causal and statistical inference without requiring the number of dependent units
to be fixed as sample size increases. We are not aware of any previous methods for inference about

2

network members in observational settings that allow the number of ties per node to increase as the
network grows. As we discuss in Section 2.3, below, this is a crucial feature of most realistic models
for social network generation. While previous methods have generally implicitly focused on one of
two possible sources of dependence among social network observations, we allow for both dependence
due to contagion or transmission of information across network ties, and dependence due to latent
similarities among nodes sharing ties. We describe estimation and inference for causal effects that are
specifically of interest in social network settings.
In Section (2) we give some background on causal inference for social network data, discussing
briefly the relationship between causal structural equation models and network edges, the types of statistical dependence likely to be found in social network data, and asymptotic growth. In Section (3.1)
we present our target of inference and the identifying assumptions that we will use in the methods that
follow. We present the efficient influence function for our target parameter under the conditional independence assumptions from van der Laan (2014). When these independence assumptions are relaxed,
this will still be an influence function for our target parameter but it may not be efficient. In Section
(3.3) we describe estimation procedures–using the (efficient) influence function from (3.1)–which will
be efficient under the stronger independence assumptions but still consistent and asymptotically normal under the weaker independence assumptions. In Section (3.4) we prove our main result, which is
the asymptotic normality of our estimator under an asymptotic regime in which the number of ties per
node grows unboundedly with n. In Section (4) we discuss estimation of rather esoteric causal effects
that are specifically of interest in social network settings. Section (5) includes simulation results, and
Section (6) concludes.

2

Background and setting

2.1

Networks and structural equation models

A network is a collection of units and information about the presence or absence of pairwise ties
between them. The presence of a tie between two units indicates that the units share some kind of
a relationship; what types of relationships are encoded by network ties depends on the context. For
example, in a social network we might define a tie to include familial relatedness, friendship, or shared
place of work. In the study of sexually transmitted diseases, such as HIV, the sexual contact network
is of great interest; here ties represent sexual partnership within a given time frame. Some types
of relationships are mutual, for example familial relatedness and shared place of work. Others, like
friendship, can go in only one direction: Tom may consider Sue to be his friend, while Sue does not
consider Tom to be her friend. Edges can be binary (present or absent), multiplex (categorical with
different levels for different types of relationships), or weighted (with continuous or ordinal weights)
according to the strength of the relationship.
Networks are most often studied through the lens of graph theory. Units are represented by nodes
and ties are represented by edges. Two nodes connected by an edge are called neighbors. A node whose
characteristics we wish to explain or model is called an ego; the ego’s neighbors are its alters. If ties
must be mutual then the graph is undirected ; if a tie can go in only one direction then the graph is
directed. For simplicity we will assume all networks are undirected in what follows, but our methods are
3

equally applicable to directed networks. In an undirected network, the degree of a node is the number
of edges it has or, equivalently, the number of alters. The graph theoretic approach to networks is
immensely powerful and useful. It has been used to study all manner of network structures, network
generating models, and diffusion processes on networks (see Newman, 2009 and references therein for
examples). However, this approach largely reduces networks to topological features, and it is not well
suited for statistical and causal inquiries into mechanisms of diffusion and influence across network
nodes.
Underlying inquiries into causal effects across network nodes is a representation of the network as a
structural equation model. Consider a network of n subjects, indexed by i, with binary undirected ties
Aij ≡ I {subjects i and j share a tie}. The matrix A with entries Aij is the adjacency matrix for the
network. Associated with each subject is a vector of random variables, Oi , including an outcome Yi ,
covariates Ci , and an exposure or treatment variable Xi , all possibly indexed by time t. In numerous
applications across the social, political, and health sciences, researchers are interested in ascertaining
the presence of and estimating causal interactions across alter-ego pairs. Is there interference, i.e.
does the treatment of subject i have a causal effect on the outcome of subject j when i and j share a
network tie? Is there peer influence, i.e. does the outcome of subject i at time t have a causal effect
on a future outcome of subject j when i and j are adjacent in the network? These inquiries can be
formalized with the help of a causal structural equation model, informed by the network.
A structural equation model is a system of equations of the form yi = fi [pai (Y ), εi ], where pai (Y ),
the set of parents of Y , is a collection of variables that are causes of Y for subject i, and εi is an error
term that may include omitted causes of Y . In general Ci and Xi will be included in pai (Y ). See
Pearl (2000) for further discussion of SEMs. When causal inference is performed on network data, the
network ties inform which variables are to be included in pai (Y ). For example, if interference might
be present, then the collection of treatment variables for i’s alters, {Xj : Aij = 1}, must be included
in the set pai (Y ). If contagion might be present then {Yj,t−k : Aij = 1} must be included in the set
pai (Yt ), where t indexes time and k is an outcome-specific lag time such that no causal effect can be
transmitted from one person to another in less than k time steps (see Ogburn and VanderWeele, 2013
for discussion of lag times).
It is important that the network be completely and accurately specified; missing ties are akin to
missing components of a multidimensional treatment vector because they result in important elements
of exposure of interest being left out of the SEM. Whenever an inquiry into causal effects is informed
by a social network, measurement error in the network is tantamount to measurement error in the
exposure of interest. The network for which data is collected must be calibrated to the causal question
of interest. If we are interested in peer effects on academic achievement among elementary school
children and think that being in the same classroom is the relationship that determines whether or
not two children affect one another’s outcomes, then being in the same classroom is the relationship
that determines whether or not a network tie exists, and a network that captures interaction during
playground sports is not informative or useful. In other words, a tie between nodes i and j represents
the possibility of a causal effect of an element of Oi on an element of Oj at a later time, and vice
versa. These issues have not been made explicit in much of the existing literature on causal inference
for network data; equating a network with the underlying SEM can help to make them precise.

4

2.2

Networks and dependence

Perhaps the greatest challenge and barrier to causal and statistical inference using observations from a
single, interconnected social network is dependence among observations. The literature on statistics for
dependent data is vast and multifaceted, but very little has been written on the dependence that arrises
when observations are sampled from a single network. Most of the literature on dependent random
variables assumes that the domain from which observations are sampled (e.g. time or geographic
space) has an underlying Euclidean topology. The principles behind asymptotic results in the Euclidean
dependence literature are simple and intuitive. They rely on a combination of stationarity assumptions,
i.e. assumptions that certain features of the data generating process do not depend on an observation’s
location in the sample domain, and assumptions that bound the nature and the amount of dependence
in the data. Most frequently these are mixing assumptions, which describe the decay of the correlation
between observations as a function of the distance between them. Sometimes the stronger assumption
of m-dependence is made, according to which two observations are independent if they are sampled
from locations that are m or more units apart. Intuitively, in order to extract an increasing amount of
information from a growing sample of dependent observations, old observations must be predictive of
new observations, which is ensured by stationarity assumptions, and the amount of independence in
the sample must grow faster than the amount of dependence, which is ensured by mixing conditions
or m-dependence.
This literature is not immediately applicable to the network setting. Roughly, this is due to the
difference between Euclidean and network topology. While it is possible to embed a network in Rd
in such a way that preserves distances, to do so is to allow d to increase as n increases. Euclidean
dependence results generally require d to be fixed, implying that, as new observations are sampled at the
boundary of a Euclidean domain, the average and maximum pairwise distance between observations
increases. Networks, on the other hand, often do not have a clear boundary to which we can add
observations in such a way that ensures growth in the sample domain. In a large sample with Euclidean
dependence, most observations will be distant from most other observations. This is not necessarily the
case in networks. The maximum distance between two nodes can be small even in very large networks,
and even if the maximum distance between two nodes is large, there may be many nodes that are close
to one another. Even under m-dependence with small m, networks exist in which most observations
remain less than m units apart as n → ∞. Therefore, mixing conditions and m-dependence do not
necessarily result in more independence than dependence in a large sample from a network. Research
indicates that social networks generally have the small-world property (sometimes referred to as the
“six degrees of separation” property), meaning that the average distance between two nodes is small
(Watts and Strogatz, 1998). Therefore distances in real-world networks may grow slowly with sample
size. Of course some types of networks, e.g. lattices, embed in Rd as n grows, but these are generally
trivial cases that are not useful for naturally occurring networks like social networks.
Dependence in networks is of two varieties, each with its own implications for inference. In the
literature on spatial and temporal dependence, dependence is often implicitly assumed to be the result
of latent traits that are more similar for observations that are close in Euclidean distance than for
distant observations. This type of dependence is likely to be present in many network contexts as well.
In networks, edges present opportunities to transmit traits or information, and this direct transmission

5

is an important additional source of dependence that depends on the underlying network structure.
Latent variable dependence will be present in data sampled from a network whenever observations
from nodes that are close to one another are more likely to share unmeasured traits than are observations from distant nodes. Homophily, or the tendency of people who share similar traits to form
network ties, is a paradigmatic example of latent variable dependence. If the outcome under study in
a social network has a genetic component, then we would expect latent variable dependence due the
fact that family members, who share latent genetic traits, are more likely to be close in social distance
than people who are unrelated. If the outcome were affected by geography or physical environment,
latent variable dependence could arise because people who live close to one another are more likely to
be friends than those who are geographically distant. Of course, these traits can create dependence
whether they are latent or observed. But if they are observed then conditioning on them renders
observations independent; therefore the methodological challenges are greater when they are latent.
Just like in the spatial and temporal dependence context, there is often little reason to think that we
could identify, let alone measure, all of these sources of dependence. In order to make any progress
towards valid inference in the presence of latent trait dependence, some structure must be assumed,
namely that the range of influence of the latent traits is primarily local in the network and that any
long-range effects are negligible. In a structural equation model, latent trait dependence would be
captured by dependence among the error terms across subjects.
Whenever one subject’s treatments, outcomes, or covariates affect other subjects’ treatments, outcomes, or covariates, observations across subjects are dependent. However, this kind of dependence,
which arrises from causal effects between subjects, has structure lacking in latent trait dependence.
Figure 1 depicts contagion in a network of three individuals. This diagram is the directed acyclic
graph representation (Pearl, 1995; Ogburn and VanderWeele, 2013) of the following structural equation model: At each time t, Yit is affected by i’s own past outcomes and those of i’s social contacts.
Individual 2 shares ties with 1 and 3 but individuals 1 and 3 are not connected. This structure implies
conditional independences: because any transmission from individual 1 to 3 must pass through 2 we
have that Y1t−2 ⊥ Y3t | Y2t−1 ; because information cannot be transmitted instantaneously we have that
Y1t−2 ⊥ Y2t−2 . If observations are observed at closely spaced time intervals then these conditional
independences can be harnessed for inference. There is no reason to think that any such conditional
independences would hold with latent variable dependence.
In this paper, we accommodate both dependence due to direct transmission and dependence due
to latent traits. We assume that both kinds of dependence are limited to dependence neighborhoods
determined by the underlying social network: each subject, or node, i can directly transmit to information, outcomes, or exposures to the nodes with which i shares a network tie, and each node i can share
latent traits with the nodes with which i shares a network tie or a mutual connection. That a subject
can only transmit to his or her immediate social contacts may be a reasonable assumption (indeed,
we can define network ties in such a way as to make this true by definition), but it is likely unrealistic
to assume that latent variable dependence only affects nodes at a distance of one or two ties, as we
assume throughout. This represents a first step towards valid statistical and causal inference under
more realistic assumptions than have been required by previous work, but future work is needed to
address more realistic–i.e. longer range–forms of latent variable dependence.

6

Figure 1: Dependence due to direct transmission

7

2.3

A note on asymptotic growth

There are many complex issues surrounding asymptotic growth of networks (e.g. Diaconis and Janson,
2007; Shalizi and Rinaldo, 2013), and a large literature on graph limits (Lovász, 2012). These issues are
largely beyond the scope of this paper, but we believe that our methods are consistent with the large
and realistic class of network-generating models.In particular, observed social networks and models
proposed for generating social networks tend to have heavy-tailed degree distributions, with most
nodes having low degree but a non-trivial proportion of nodes having high degree, with the maximum
degree dependent on the size of the network. Some researchers speculate that the heavy right tails of
social network degree distributions tend to approximately follow power laws: P r(degree = k) ∼ k −α
for 2 < α < 3 (Barabási and Albert, 1999; Lovász, 2012; Newman and Park, 2003), in which case
P r(degree > K) = O(K 1−α ) for any fixed K. Even if degree distributions depart from power law
distributions (Clauset et al., 2009) they are frequently incompatible with the assumption of bounded
degrees, which has been used in previous methods for inference about observations sampled from a
single social network. Our new methods are not able to accommodate the most highly connected nodes
from a power law degree sequence, but they can nevertheless be used to perform inference about the
other nodes in a network that has a power law degree distribution (see Section 4.4).
Our theoretical results require an asymptotic regime in which the number of nodes in the network,
n, goes to infinity. Formalizing asymptotic growth of network-generating models is beyond the scope
of this paper; we take for granted sequence of networks with increasing n such that key features of the
network–whatever features inform the causal effects of interest–are preserved. The structural equation
model that specifies the distributions of covariates, treatment, and outcome does not change as n
increases.

3

Methods

In this section we describe estimation of and inference about the causal effect of a treatment or
exposure, X, which can depend upon covariates but is otherwise independent of outcomes and network
topology. This case includes randomized and non-randomized exposures that may or may not be
subject to interference, specifically direct interference (Ogburn and VanderWeele, 2013), where the
effect of one individual’s exposure on another’s outcome is unmediated by the first individual’s outcome.
(In Section 4 we discuss causal effects that correspond to interference due to contagion, where the effect
of one individual’s exposure on another’s outcome is mediated by the first individual’s outcome, and
allocational interference, where treatment changes network ties.) The approach we describe below is
different from traditional approaches to interference in that they are justified when partial interference
does not hold. As far as we are aware, this is the first approach to interference that references an
asymptotic regime in which the number of ties for a given individual may grow with sample size. The
estimating procedure that we describe in this section was first proposed by van der Laan (2014), but
we generalize the results to a broader class of causal effects and to more general and pervasive forms of
dependence among observations. The conditions under which the resulting estimators are consistent
and asymptotically Normally distributed are different and weaker here than those in van der Laan
(2014).

8

Table 1: Properties of marginal estimands and of estimands conditional on C
Properties that we have demonstrated for the two classes of estimands
nonparametrically identified with or without latent variable (LV) dependence
estimator is CAN with or without LV dependence
efficient estimator is available with LV dependence
efficient estimator is available without LV dependence
consistent and computationally tractable variance estimation with LV dependence
consistent and computationally tractable variance estimation without LV dependence

Estimand class
Marginal Conditional
yes
yes
yes
yes
no
no
yes
yes
no
yes
yes
yes

For the remainder of Section 3, we describe CAN estimators of causal effects under two different sets
of assumptions. One set of assumptions allows dependence due to direct transmission but not latent
variable dependence, as in van der Laan (2014); under this set of assumptions our estimators inherit
the efficiency properties from van der Laan (2014). The other set of assumptions allows dependence
due to direct transmission and latent variable dependence; under this set of assumptions our estimators
are CAN but may not be efficient. Our main result is the proof of asymptotic normality under an
asymptotic regime in which the number of ties for a given individual may grow with sample size in
Section 3.4.
In Section 3.5 we describe statistical inference for the estimators introduced in Section 3. We
consider two different classes of estimators: estimators that marginalize over baseline covariates and
estimators that condition on baseline covariates. In some cases, variance estimation is facilitated by
conditioning on covariates. Under the assumptions encoded in the structural equation model in Section
3.1, the conditional estimator is in fact consistent for the marginal estimand. However, conditional
estimators have smaller variance and inference about the conditional estimand cannot be interpreted as
inference about the marginal estimand. All of our estimands and estimators condition on the observed
network as given by the adjacency matrix A. Table 1 summarizes the relationships among the two
sets of assumptions (with and without latent variable dependence) and the two classes of estimators
(marginal over C and conditional on C) according to their properties and according to the limitations
of our proposed methods.
In 4 we describe more elaborate interventions and causal questions that can be addressed by our
methods and that may be of interest specifically in the context of social network data. However,
we focus throughout on single treatment. Longitudinal interventions are also possible but we leave
the details for future work. In addition, we state our results under the assumption that all variables
take values on discrete sets. Analogous results are valid for other types of random variables: it
is straightforward to extend our notation and results to continuous covariates and outcomes, but
continuous treatments are more complicated (see van der Laan, 2014).

3.1

Structural equation model

Recall that Ci denotes a vector of covariates, Xi is a treatment variable, and Yi is the outcome for
Pn
subjects i ∈ {1, . . . , n}. In addition, let Ki = j=1 Aij , that is, Ki is degree of node i or the number
of individuals sharing a tie with individual i. The degree of subject i and the degrees of i’s alters may
be included in the covariate vector Ci . We define Y = (Y1 , ..., Yn ) and C and X analogously. We
9

use a structural equation model to define the causal effects of interest, as in Section 2, but note that
analogous definitions may be achieved within the potential outcome framework (van der Laan, 2014).
We assume that the data are generated by sequentially evaluating the following set of equations:
Ci = fC [εCi ]

i = 1, . . . , n

Xi = fX [{Cj : Aij = 1} , εXi ]

i = 1, . . . , n

Yi = fY [{Xj : Aij = 1} , {Cj : Aij = 1} , εYi ]

i = 1, . . . , n,

(1)

where fC , fX , and fY are unknown and unspecified functions and εi = (εCi , εXi , εYi ) is a vector of
exogenous, unobserved errors for individual i. This set of equations encodes the assumptions of direct
interference in an observational or randomized setting: the treatment X may or may not depend
on C, depending on the specification of fX , and the outcome Y depends on X and C but not on
previous outcomes or on A, except possibly through C. Time ordering is a fundamental component
of a structural causal model. For example, we assume that C is first drawn for all units, so that, in
addition to Ci , the other components of the vector C may affect the value of Xi .
In addition, nonparametric identification of causal effects requires the following assumptions on the
error terms from the SEM:

(εX1 , ..., εXn ) ⊥ (εY1 , ..., εYn ) | C,

(A1)

εX1 , ..., εXn | C are identically distributed and εY1 , ..., εYn | C, X are identically distributed, (A2a)
εXi ⊥ εXj | C and εYi ⊥ εYj | C, X for i, j s.t. Aij = 0 and ∃!k with Aik = Akj = 1

(A2b)

εCi , i = 1, ..., n, are identically distributed, and

(A3a)

εCi ⊥ εCj for i, j s.t. Aij = 0 and ∃!k with Aik = Akj = 1.

(A3b)

Assumption (A1) is a no unmeasured confounding assumption; it ensures that C suffices to control for
confounding of the effect of X on Y. It implies that any latent variable dependence affects X and Y
separately; in general a latent variable that affected X and Y jointly would constitute a violation of
this assumption. Assumptions (A2b) and (A3b) ensure that any unmeasured sources of dependence–
i.e. latent trait dependence–only affect pairs of observations up to a distance of two network ties–that
is, friends or friends-of-friends. Assumption (A3) can be omitted if attention is restricted to causal
effects conditional on C.
Although our main result, given in Theorem 1 below, holds for inference in the SEM defined by
assumptions (A1)–(A3b), some asymptotic properties are guaranteed only when stronger versions of
assumptions (A2b) and (A3b) hold. We therefore introduce alternative assumptions
εX1 , ..., εXn | C are i.i.d. and εY1 , ..., εYn | C, X are i.i.d., and

(A4)

εCi , i = 1, ..., n, are i.i.d.

(A5)

Note that, although the variance-covariance structure of the SEM given in (1) is affected by the
dependence allowed in (A2b) and (A3b), the mean structure is unaltered by the choice of assumptions

10

(A2) and (A3) or (A4) and (A5). Therefore, any estimator that is unbiased under (A4) and (A5) will
remain unbiased when these are relaxed to (A2) and (A3). In Section 3.2 we discuss nonparametric
identification of causal parameters, which is agnostic to the choice of the weaker or stronger independence assumptions. In Section 3.3 we derive estimators under assumptions (A1), (A4), and (A5)–that
is, in the presence of dependence due to direct transmission but not latent variable dependence. We
use the stronger assumptions because the resulting model is amenable to familiar tools for deriving
semiparametric estimators. In Section (3.4) we prove that the estimators derived under assumptions
(A1), (A4), and (A5) are CAN under the weaker set of assumptions (A1)–(A3b). In Section 3.5 we
discuss inference under each of the two sets of assumptions.

3.2

Definition and nonparametric identification of causal effects

In principle it is possible to perform statistical inference in the model defined by assumptions (A1)(A3b) or by assumptions (A1), (A4), and (A5). However, in practice, with finite data and finite computing resources, we may need to make dimension-reducing assumptions on the forms of fX and fY . This
is done by considering summary functions sX and sY and random variables Wi = sX,i ({Cj : Aij = 1})
and Vi = sY,i ({Cj : Aij = 1} , {Xj : Aij = 1}) such that the model may be written as
Ci = fC [εCi ]

i = 1, . . . , n

Xi = fX [Wi , εXi ]

i = 1, . . . , n

Yi = fY [Vi , εYi ]

i = 1, . . . , n.

For example, sX,i ({Cj : Aij = 1}) =



Ci ,

P

j:Aij =1

Cj



implies that the treatment status of unit i

only depends on i’s own covariate value and on the sum of the covariate
units sharing
 P values of theP
 a
tie with i. Analogously, sY,i ({Cj : Aij = 1} , {Xj : Aij = 1}) = Ci , j:Aij =1 Cj , Xi , j:Aij =1 Xj is
an example of a summary function for fY . For convenience we use the notation sX,i (C) and sY,i (C, X)
below; however, this notation should not undermine the important fact that Wi can only depend on the
subset and {Cj : Aij = 1} and Vi can only depend on the subsets {Cj : Aij = 1} and {Xj : Aij = 1}
of C and X, as these are the only components of C and X that are parents of X and Y, respectively,
in the causal model defined by the network-as-structural-causal-model. For notational convenience,
in what follows we augment the observed data random vector with Vi and Wi , recognizing that these
are deterministic functionals of Ci and Xi , defined by sY,i and sX,i , and are therefore technically
redundant: Oi = (Ci , Wi , Xi , Vi , Yi ).
Potential or counterfactual outcomes can be defined in terms of hypothetical interventions to the
underlying structural causal model. For example, a hypothetical intervention that deterministically
sets Xi to a user-given value x∗i for i = 1, ..., n is given by
Ci = fC [εCi ]

i = 1, . . . , n

x∗i

i = 1, . . . , n

Xi =

Yi (x∗ ) = fY [Vi (x∗ ), εYi ]

i = 1, . . . , n,

where x∗ = (x∗1 , . . . , x∗n ). Here Yi (x∗ ) denotes the potential outcome of individual i in a hypothetical
11

world in which P (X = x∗ ) = 1. Analogously, Vi (x∗ ) = sY,i (C, x∗ ) is a counterfactual random variable
in a hypothetical world in which P (X = x∗ ) = 1. Note that, although Vi (x∗ ) is counterfactual, its
value is determined by the observed realization of C and by the user-specified value x∗ , and it is
therefore known. In order to streamline notation as we describe increasingly complex interventions,
we denote the counterfactual variables Vi (x∗ ) and Yi (x∗ ) by Vi∗ and Yi∗ , respectively. The causal
parameter of interest throughout is the expected average potential outcome in this same hypothetical
 
Pn
world, i.e. E Ȳ ∗ , where Ȳ ∗ = n1 i=1 Yi∗ .
We are now ready define notation that we will use throughout the remainder of the paper for
functionals of the distribution of O. Let pC (c) = P (C = c), g(x|w) = P (X = x | W = w), gi (x|w) =
P (Xi = x|Wi = w), pY (y|v) = P (Y = y|V = v), and pY,i (y|v) = P (Yi = y|Vi = v). Define the two
marginal distributions hi (v) = P (Vi = v) and hi,x∗ (v) = P (Vi∗ = v), noting that both hi and hi,x∗
P
are determined by g and pC and are therefore observed data quantities. Finally, m(v) = y y pY (y|v)
is the conditional expectation of Y given V = v.
 
In addition to assumptions (A1)-(A3b) or (A1), (A4), and (A5), identification of E Ȳ ∗ requires
the positivity assumption that
P (V = v|C = c) > 0 for all v in the range of Vi∗ and for all c such that P (C = c) > 0.

(A6)

This assumption states that, within levels of C, the values of V determined by the hypothetical
intervention x∗ have positive probability under the observed-data-generating distribution. Now the
 
causal parameter E Ȳ ∗ is identified by
n

ψ=

n

1 XX
1X
E [m(Vi∗ )] =
m(v)hi,x∗ (v).
n i=1
n i=1 v

(2)

This identification result is equivalent to
n

ψ=

1 XX
m(sY,i (c, x∗ ))pC (c).
n i=1 c



From (3), it is clear that the conditional causal parameter E Ȳ ∗ | C = c is identified by

3.3

(3)
1
n

Pn

i=1

m(sY,i (c, x∗ )).

Estimation

 
Estimation of and inference about E Ȳ ∗ requires a statistical model M for the distribution of the
observed data P (O). That is, M is a collection of distributions over O of which one element is
the true data-generating distribution. Our target of inference is a pathwise differentiable mapping
Ψ : M → R such that ψ is Ψ(P ), the mapping evaluated at the true data-generating distribution.
Under assumptions (A1), (A4), and (A5) the probability distribution of the observed data may be
factorized as
P (O = o) = P (C = c) g(x|w)pY (y|v),

(4)

suggesting that M requires three components: a model for pC , a model for g, and a model for P (Y |V ).
Furthermore, the identification results in (2) and (3) indicate that ψ depends on P (Y |V ) only through

12

m. The empirical distribution p̂C can be used throughout to nonparametrically estimate pC , but,
when C is high-dimensional, g and m cannot be non-parametrically estimated at rates of convergence
that are fast enough to satisfy the regularity conditions of Theorem 1 (see Appendix). Therefore, in
order to define the parameter mapping we require a statistical model M = Mg × Mm , where Mg is
a collection of conditional distributions for X given W such that the true conditional distribution is a
member, and Mm is a collection of expectations of Y relative to conditional distributions of Y given
V such that the true conditional expectation of Y given V is a member. Estimation of ψ is based on
the efficient influence function for the parameter mapping Ψ : M → R.
Under assumptions (A1), (A4), and (A5), the efficient influence function, D, evaluated at a fixed
value o of O, is given by
D(o) =
where h̄(vi ) =

1
n

Pn

n
n
n
X
1X
1 X h̄x∗ (vi )
E [m (Vi∗ ) | Cj = cj ] − ψ +
{yi − m (vi )} ,
n i=1
n i=1 h̄(vi )
j=1

j=1

hj (vi ), h̄x∗ (vi ) =

1
n

Pn

j=1

(5)

hj,x∗ (vi ), vi = sY,i (c, x), and Vi∗ = sY,i (C, x∗ ). The

influence function has expected value equal to 0 at the true ψ; this fact can be used to generate
unbiased estimating equations for ψ. Estimating equations based on the efficient influence function
are doubly robust: the right hand side of Equation (5) has expected value equal to 0 if m(·) is replaced
with an arbitrary functional of V or if g(·) is replaced with an arbitrary functional of W , as long as
one of the two remains correctly specified. (Recall that g(·), along with pC , determines h̄x∗ (vi ) and
h̄(vi ).) This implies that an estimating equation based on Equation (5) will be unbiased for ψ if either
model Mm for m(·) or model Mg for g(·) is correctly specified, i.e. contains the truth, even if one is
not. This influence function is efficient in that, when m(·) is correctly specified, it has the smallest
variance among all influence functions in model Mg .
The efficient influence function in a model that does not make any distributional assumptions about
C, that is under assumptions (A1) and (A4) only, is given in equation (6) below.

n 
h̄x∗ (vi )
1X
∗
E [m (Vi ) | C = c] − ψ +
{yi − m (vi )} .
D (o) =
n i=1
h̄(vi )
0

(6)

We use this influence function in the estimating procedure below. This is also the influence function
used to derive estimators conditional on C, in which case the first two terms cancel out; we will denote
the conditional influence function wit Dc (o).
In social network settings C is likely to be high-dimensional – recall that Ci includes all relevant
covariates for individual i and for i’s social contacts, in addition to a list of network ties. Even though
we assume that m(·) and g(·) depend on C only through the dimension-reducing functions sY (·)
and sX (·), if m(·) and g(·) are estimated nonparametrically care should be taken to ensure rates of
convergence satisfying the regularity conditions of Theorem 1. When m(·) and g(·) are estimated with
working models, the doubly robust property is advantageous, affording two independent opportunities
for valid inference about ψ. If models for both m(·) and g(·) are correctly specified, then the solution
to estimating equations based on the efficient influence function are asymptotically efficient in the
semiparametric model defined by the identifying assumptions and by the correctly specified model for
g(·). See Van der Laan and Robins (2003); Tsiatis (2007) for a review of semiparametric inference in

13

i.i.d. settings.
A traditional estimating equation approach to inference about ψ has two disadvantages in finite
samples. First, it does not result in substitution estimator, leading to the possibility that finite-sample
estimates will fall outside of the parameter space when the parameter space for ψ is bounded. Second,
it is sensitive to extreme values of the term h̄i,x∗ (v)/h̄(vi ) in the righthand side of Equation (5).
For these reasons, we propose a targeted maximum loss-based estimator (TMLE) of ψ; TMLEs are
substitution estimators and are not as sensitive to the near violations of the positivity assumption that
can occur in finite samples and result in extreme values of h̄x∗ (vi )/h̄(vi ).
Targeted maximum likelihood estimation is a general template for estimation of smooth parameters
in semi- and nonparametric models. The estimation algorithm is constructed to solve the efficient
influence function estimating equation, thereby yielding, under regularity conditions, asymptotically
linear estimators with the same semiparametric efficiency property as the estimating equation approach
described above. In our setting, a TMLE is constructed using three elements: (i) a valid loss function
L for the outcome regression model m, (ii) initial working estimators m̂ of m and and ĝ of g, and
(iii) a parametric submodel m of M, the score of which corresponds to a particular component of
the score based on the efficient influence function D(o) and such that m0 = m(·). The TMLE is then
defined by an iterative procedure that, at each step, estimates  by minimizing the empirical risk of
the loss function L at m . An updated estimate is then computed as m̂ˆ, and the process is repeated
until convergence. The TMLE is the estimator obtained in the final step of the iteration. The result
of the previous iterative procedure is that, at the final step, the efficient influence function estimating
equation is solved. For more details about targeted maximum likelihood estimation, see Van der Laan
and Rose (2011). Below we describe a TMLE for ψ based on (6) that requires only one iteration for
convergence. We use influence function D0 (o) to derive the TMLE, instead of D(o), because it is
computationally more tractable and because the choice of influence function does not matter for the
conditional parameter that we are interested in when latent variable dependence is present.
In order to use targeted maximum likelihood estimation to estimate ψ, the parameter identifying
 
E Ȳ ∗ , initial estimators m̂ and ĝ of m and g may be found through standard maximum likelihood
or loss-based estimation methods. For example, if Yi is binary, or if it is continuous and bounded and
rescaled to fall between 0 and 1, then the negative log likelihood function
n
X


log m(Vi )Yi (1 − m(Vi ))1−Yi

i=1

is a valid loss function for m. Assuming a parametric model mθ (·) for m(·), a maximum likelihood
estimator is given by mθ̂ , where
θ̂ = arg max
θ

n
X


log mθ (Vi )Yi (1 − mθ (Vi ))1−Yi

i=1

may be estimated by running standard binary regression software of the outcomes Yi on Vi . An initial
Pn
i=1 log g(Xi |Wi ).

estimator ĝ may be found analogously by optimizing the log likelihood function

Alternatively, these loss functions may be used to implement other loss-based estimation techniques
such as super learning (van der Laan Mark et al., 2007), which selects the weighted average of a
14

collection (“library”) of estimators that minimizes the cross-validated loss-based risk; the library may
include any machine learning algorithms and logistic regression estimators. The empirical distribution
p̂C is used to estimate pC .
Estimation of h̄ and h̄x can be carried out by substituting ĝ and p̂C for g and pC in the expressions

h̄(v) =
h̄x∗ (v) =

1 XX
I(sY,j (c, x) = v)g(x|w)pC (c) and
n j c,x
1 XX
I(sY,j (c, x∗ ) = v)pC (c),
n j c

ˆ and h̄
ˆ ∗ the corresponding
where I(E) is the indicator function of the event E. We denote by h̄
x
estimates of h̄ and h̄x∗ .
An alternative and computationally more feasible approach is to estimate h̄(v) by directly optimizPn
ing the log likelihood function i=1 log h̄(Vi |Wi ), as if the pooled sample (Vi , Wi ) were i.i.d. It can
be shown that this results in a valid loss function for h̄, even for dependent observations (Vi , Wi ), for
i = 1, . . . , n. Similarly, one can construct a direct estimator of h̄x∗ , by first creating a sample (Vi∗ , Wi )
Pn
and then directly optimizing the log likelihood function i=1 log h̄x∗ (Vi∗ |Wi ), as if the pooled sample
(Vi∗ , Wi ) were i.i.d.
Now the targeted minimum loss based estimator of ψ is computed as follows:
1. Define the auxiliary weights Hi as the ratio of estimated densities of V ∗ and V evaluated at the
observed value Vi . Compute the auxiliary weights as
Hi =

ˆ ∗ (V )
h̄
x
i
.
ˆ
h̄(V )
i

2. Compute initial predicted outcome values Ŷi = m̂(Vi ) and predicted potential outcome values
Ŷi∗ = m̂(Vi∗ ) evaluated at the counterfactual value Vi∗ = sY,i (C, x∗ ).
3. Construct a TMLE model update m̂ˆ of m̂ by running a weighted intercept-only logistic regression
model with weights Hi defined in step (1), Yi as the outcome and including Ŷi as an offset. That
is, define ˆ as the estimate of the intercept parameter  from the following weighted logistic
regression model
logitm̂ (v) = logitm̂(v) + ,
where logit(x) = log



x
1−x



.
∗

4. Compute updated predicted potential outcomes Ỹi as the fitted values of the regression from
step 3, evaluated at v ∗ rather than v (that is, at Ŷi∗ instead of Ŷi ):
∗

Ỹi = expit{logitŶi∗ + ˆ},
where expit(x) =

1
1+e−x ,

i.e., the inverse of the logit function.

15

5. Compute the TMLE ψ̂ as
n

ψ̂ =

1X ∗
Ỹi .
n i=1

The TMLE inherits the double robustness property of the estimating equation estimator we described
above: it will be consistent for ψ if either the working model ĝ for g or the working model m̂ for m
is correctly specified. This resulting estimator remains CAN for ψ under assumptions (A2) and (A3)
instead of (A4) and (A5), and the same procedure can be used to estimate the parameter conditional
on C.

3.4

Asymptotic normality

van der Laan (2014) proved the asymptotic normality of the estimators described above under assumptions that the degree Ki is uniformly bounded above as n → ∞ by a constant K. In order to
accommodate more realistic models of asymptotic growth in the network context, we consider a new
asymptotic regime in which Ki may grow as n → ∞.
2
Theorem 1: Let Kmax,n = maxi {Ki } for a fixed network with n nodes. Suppose that Kmax,n
/n →

0 as n → ∞. If at least one of the two models for m(·) and g(·) is correctly specified and under
independence assumptions (A1) through (A3b), positivity assumption (A6), and regularity conditions
(see Appendix),

p 
d
Cn ψ̂ − ψ −→ N (0, σ 2 ),
2
n/Kmax,n
≤ Cn ≤ n. The asymptotic variance of ψ̂, σ 2 , is given by the variance of the influence curve

of the estimator.
The proof of Theorem 1 is in the Appendix. In Section 4.4, below, we discuss settings in which the
conditions for this theorem fail to hold, and ways to recover valid inference for conditional estimands
in some of these settings. Broadly, the proof has two parts: first, to show that the second order terms
√
in the expansion of ψ̂ − ψ are stochastically less than 1/ Cn , and second, to show that the first order
√
terms converge to a normal distribution when scaled by a factor of order Cn . The proof that the
√
second order terms are stochastically less than 1/ Cn is an extension of the empirical process theory
of Van Der Vaart and Wellner (1996) and follows the same format as the proof in van der Laan (2014).
For the proof that the first order terms converge to a normal distribution, we rely on Stein’s method of
central limit theorem proof (Stein, 1972). Stein’s method allows us to derive a bound on the distance
between our first order term (properly scaled) and a standard normal distribution; this bound depends
on the degree distribution K1 , ..., Kn . We show that this bound converges to 0 as n → ∞ under
2
regularity conditions and our running assumption that Kmax,n
= o(n).

When all nodes have the same
q number of ties, i.e. Ki = Kmax,n for all i, then the rate of
√
2
convergence will be given Cn = n/Kmax,n
. When Kmax,n is bounded above as n → ∞, then the
√
rate of convergence will be n. This is the setting covered in van der Laan (2014). When Kq
max,n → ∞
2
but some nodes have fewer than Kmax,n ties, the exact rate of convergence is between n/Kmax,n
√
and n (inclusive) but is difficult or impossible to determine analytically, as it may depend intricately

on the topology of the network. The inferential procedures that we describe below do not require
16

knowledge of the rate of convergence.

3.5

Inference

√
A 95% confidence interval for ψ is given by ψ̂n ± 1.96σ/ Cn . In practice neither σ nor Cn are likely to
be known, but available variance estimation methods estimate the variance of ψ̂n directly, incorporating
the rate of convergence without requiring it to be known a priori.
In i.i.d settings, the variance of doubly robust estimators can be bootstrapped or, under the assumption that both nuisance models are correctly specified, estimated with the plug-in estimator of
the square of the estimator’s influence function. When the nuisance parameters are estimated with
parametric models, standard sandwich variance estimators are appropriate whether one or both of
the nuisance models are correctly specified. In the present, non-i.i.d. setting, variance estimation is
considerably more challenging. In principle, the variance of ψ̂ can be estimated using the empirical
average of the square of the influence function, substituting ψ̂ for ψ and the fitted values from the
working models ĝ and m̂ for g and m. Although this variance may be anticonservative if one, but
not both, of the working models ĝ and m̂ is correctly specified, using flexible or non-parametric specifications for these models increases opportunities to estimate both consistently. However, unlike in
i.i.d. settings, the expectation of the square of the empirical version of the influence function given in
Equation (5) does not reduce to the sum of squared influence terms for each observation. Instead, it
includes double sums for all pairs of observations that are not marginally independent of one another.
These terms capture covariances between dependent observations; these extra covariance terms reflect
a larger variance and a slower rate of convergence due to dependence across observations.
When dependence is due to direct transmission, that is, under assumptions (A1), (A4), and (A5),
two alternative variance estimation procedures are available. One option is to estimate the variance
of the influence function D0 (o) given by Equation (6). Our TMLE is based on D0 (o), but because this
is the efficient IF in a model that makes fewer assumptions than (A1), (A4), and (A5), it has larger
variance than D(o) and provides a valid (asymptotically conservative) variance estimate even when
estimation is based on D(o). For consistent and computationally feasible estimators for the variance
of D0 (o) see Sofrygin and van der Laan (2015).
An alternative approach to estimate the variance of ψ̂ under assumptions (A1), (A4), and (A5) is
to employ the following version of a parametric bootstrap. This approach might offer improvements
in finite-sample performance over the previously described approach. Briefly, our proposed procedure
proceeds by iteratively sampling n observations from the existing fit of the likelihood, fitting the
univariate least favorable submodel parameter , and computing the bootstrapped TMLE. After enough
iterations one can obtain the Monte Carlo variance estimator of ψ̂ by evaluating the empirical variance
of the bootstrap TMLEs. In more detail, for each of M bootstrap iterations, indexed as b = 1, . . . , M ,
first n covariates Cb = (C1b , . . . , Cnb ) are sampled with replacement (assuming C are i.i.d.), then the
existing model fit ĝ is applied to sampling of n exposures Xb = (X1b , . . . , Xnb ), followed by a sample
of n outcomes Yb = (Y1b , . . . , Ynb ) based on the existing outcome model fit m̂. Note that we are also
assuming that the corresponding bootstrap random summaries Wib and Vib , for i = 1, . . . , n, were
constructed by applying the summary functions sX and sY to Cb and (Cb , Xb ), respectively. This
bootstrap sample is then applied to obtain the predicted values from the existing auxiliary covariate

17

ˆ ∗ /h̄
ˆ )(V b ), for i = 1, . . . , n, followed by a bootstrap-based fitting of , and finally, evaluation of
fit (h̄
x
i
bootstrap TMLE. Note that the TMLE model update is the only model fitting step needed at each
iteration of the bootstrap, which significantly lowers the computational burden of this procedure. The
variance estimate is then obtained by taking the empirical variance of bootstrap TMLE samples ψ̂ b . In
the present setting, due to dependence across observations, one must be judicious with applications of
the bootstrap. For example, the parametric bootstrap procedure described above requires conditional
independence of Xi given Wi and Yi given Vi , along with the consistent modeling of the corresponding
factors of the likelihood. It may seem natural to sample Vi directly from its corresponding auxiliary
model fit, but this is likely to result in an anti-conservative variance estimates, since the conditional
independence of Vi is unlikely to hold by virtue of its construction as a summary measure of the
network.
When latent variable dependence is present, that is under assumptions (A1) through (A3), consistent and computationally feasible variance estimation procedures are not currently available for either
D0 (o) or D(o), because existing methods require bootstrapping some of the observed data. Without
latent variable dependence we can take advantage of marginal and conditional independences to employ i.i.d. or parametric bootstrap methods, but latent variable dependence requires new methods
for dependent data bootstrap. For this reason, we instead estimate the conditional parameter with
influence function DC (o). A simple plug-in estimator is available for the variance of this influence
function (see the Appendix and van der Laan, 2014).

4

Extensions

In this section we extend the estimation procedure to two causal effects of great interest in the context
of social networks: to social contagion, or peer effects; to stochastic interventions; and to interventions
on the network topology itself, i.e. interventions on A = [Aij : i, j ∈ {1, . . . , n}] where, as above,
Aij ≡ I {subjects i and j share a tie}. First we introduce dynamic and stochastic interventions.

4.1

Dynamic and stochastic interventions

Interventions that assign treatment as a user-specified function dX (·) of C correspond to substituting
dX,i (C) for x∗i in the intervention model, definitions, and estimating procedure above. This is an
example of a “dynamic” intervention: treatment is deterministically specified conditional on covariates
but is but allowed to depend (“dynamically”) on covariates.
We can also identify the effects of interventions that replace fX with a new, user-specified function.
This is an example of a stochastic intervention: the intervention changes the distribution of X but
does not eliminate the stochasticity introduced by εX ; it is represented by an intervention SEM that
replaces the equation for Xi with Xi∗ = rX [Wi∗ , εXi ] for a user specified function rX . For discussions of
stochastic interventions and their identifying assumptions in i.i.d. settings, see Muñoz and van der Laan
(2012); Haneuse and Rotnitzky (2013); Young et al. (2014). In the social network setting, stochastic
interventions that change the dependence of Xi on C and of and Yi on C and X are of particular
interest. For example, consider data generated by an SEM in which fX depends on C only through
P
Wi = |A1i | j:Aij =1 Cj , i.e. the mean of C among the set of alters of i. We might be interested in
18

the mean counterfactual outcome under a stochastic intervention that forces fX to depend instead on
Wi∗ = maxj:Aij =1 {Cj }, i.e. the maximum value C among the alters of i. This particular stochastic
intervention modifies fX only through W ; it is represented by an intervention SEM that replaces the
equation for Xi with Xi∗ = fX [Wi∗ , εXi ]. This is not a deterministic intervention, rather it assigns X
according to a stochastic distribution: for each x in the support of X, Xi is set by the intervention to


x with probability P X = x|W = maxj:Aij =1 {Cj } .
We formally define the class of stochastic interventions that alter the dependence of Xi on C and
of and Yi on (C, X), discuss identifying assumptions and estimation procedures, and then describe
some such interventions of particular interest. Let s∗X,i (·) and s∗Y,i (·, ·) be user-specified functionals.
They are denoted by an asterisk because they index hypothetical interventions rather than realized
data-generating mechanisms. Let Wi∗ = s∗X,i (C) and Vi∗ = s∗Y,i (C, X∗ ). We are concerned with the
class of stochastic interventions given by
Ci = fC [εCi ]

i = 1, . . . , n

Xi∗ = fX [Wi∗ , εXi ]

i = 1, . . . , n

Yi∗

i = 1, . . . , n.

=

fY [Vi∗ , εYi ]

(7)

This can be interpreted as an intervention where, for each x∗ in the support of X and for i = 1, ..., n,


Xi is set to x∗ with probability P X = x∗ |W = s∗X,i (C) and Vi is set to s∗Y,i (C, x∗ ) deterministically
for each possible realization x∗ . Because Y depends on X only through V , this is equivalent to an



intervention that sets Vi to v with probability P X ∈ x∗ : s∗Y,i (C, x∗ ) = v | W = s∗X (C) , where

s∗X (C) = s∗X,1 (C), ..., s∗X,n (C) .
This intervention is identified under the same assumptions as the deterministic interventions described above, with the exception of a positivity assumption that is a slight modification of (A6).
Define X ∗ = {x∗ : P [X = x∗ |W = s∗X (C)] > 0} to be the set of treatment vectors x∗ that have positive probability under the stochastic intervention defined by (7). We assume that

minv∈V ∗ P (V = v|C = c) > 0 for V ∗ = s∗Y,i (C, x∗ ) : x∗ ∈ X ∗ and for all c in the support of C. (8)
Note that, in order for this positivity assumption to hold, the supports of s∗X (·) and s∗Y (·, ·) must be
of the same dimensions as the supports of sX (·) and sY (·, ·), respectively.
The causal parameter of interest is the expected average potential outcome under this hypothetical
 


 
intervention, E Ȳ ∗ . Define h∗i (v) = P [Vi∗ = v] = P s∗Y,i (C, X∗ ) = v . Then E Ȳ ∗ is identified by
n

ψ=


1 XX 
E Yi |s∗Y,i (c, x) P [X = x|W = s∗X (c)] pC (c)
n i=1 c,x
n

=

n

1X
1 XX
E [m(Vi∗ )] =
m(v)h∗i (v).
n i=1
n i=1 v

19

An influence function for ψ, evaluated at a fixed value of the observed data, o, is given by
D(o) =

where h̄∗ (vi ) =

1
n

n
n
n
X
1X
1 X h̄∗ (vi )
{yi − m (vi )} ,
E [m (Vi∗ ) | Cj = cj ] − ψ +
n i=1
n i=1 h̄(vi )
j=1

Pn

j=1
∗

h∗j (vi ). (This is the efficient influence function under assumptions (A4) and

(A5).) Estimation of h̄ is carried out by substituting ĝ and p̂C for g and pC in the expression
h̄∗ (v) =


1 XX
I s∗Y,i (c, x) = v g (x|s∗X (c)) pC (c).
n j c,x

ˆ∗
Since p̂C is an empirical distribution that puts mass one on the observed value c, the estimator h̄
reduces to
ˆ ∗ (v) = 1
h̄
n

n X
X


I s∗Y,i (x, C) = v ĝ(x|s∗X (C)).

j=1 x

ˆ and h̄
ˆ ∗ the corresponding estimates of h̄ and h̄∗ .
We denote by h̄
Now the targeted minimum loss based estimator of ψ is computed according to the steps outlined
in Section 3, but with V ∗ and Y ∗ defined as immediately above.
A special case of this class of stochastic interventions intervenes only on sX , like the example
discussed above in which the intervention forces fX to depend on Wi∗ = maxj:Aij =1 {Cj } but does not
alter the functional form of sY . E[Ȳ ∗ ] under this type of intervention is identified by
n

ψ=

1 XX
E [Yi |C = c, X = x] P [X = x|W = s∗X (c)] pC (c)
n i=1 c,x
n

=

n

1 XX
1X
E [m(Vi∗ )] =
m(v)h∗i (v).
n i=1
n i=1 v

With Vi∗ defined as sY,i (C, X∗ ), estimation of this class of intervention proceeds as immediately above.
The fact that X∗ is random does not affect the estimation algorithm.

4.2

Peer effects

Define Yi0 to be the outcome variable measured at a time previous to the primary outcome measurement
Yi . Peer effects are the class of causal effects of Yj0 on Yi for Aij = 1: the effects of individuals’ outcomes
on the subsequent outcome of their alters. We can operationalize peer effects as the effects of dynamic
interventions where the counterfactual exposure for subject i is given by a user-specified function
dX (·) of {Yj0 : Aij = 1}. In order to maintain the identifying assumptions A2b and A3b, the time
elapsed between Y 0 and Y must permit transmission only between nodes and their immediate alters.
Otherwise, if between Y 0 and Y the outcome could have spread contagiously more broadly, there will
be more dependence present than our methods can account for, and also possible confounding of the
effect of Yi0 on Yj for Aij = 1 due to mutual friends.

20

4.3

Interventions on the network topology

An intervention on the network, i.e. an intervention that adds, removes, or relocates ties in the
network, is a special case of a joint intervention on sX (·) and sY (·). To see this, note that the network
topology, codified by the adjacency matrix A, enters the data-generating structural equation model (1)
only through sX (·) and sY (·); therefore we can represent any modification to A via the corresponding
modification to sX (·) and sY (·). Consider an intervention that replaces the observed adjacency matrix
A with a user-specified adjacency matrix A∗ . This is an example of the stochastic interventions


∗
described on page 19, with s∗X,i (C) replaced by sA
Cj : A∗ij = 1 and s∗Y,i (C, X∗ ) by
X,i (C) ≡ sX,i
 ∗


∗
∗
sA
Xj : A∗ij = 1 , Cj : A∗ij = 1 . The intervention SEM differs from the dataY,i (C, X ) ≡ sY,i
generating SEM only in that Xi depends on the covariate values for the individuals with whom i
shares ties in the intervention adjacency matrix A∗ and Yi depends on the counterfactual treatments
and observed covariate values for those same individuals.
Interventions on summary features of the adjacency matrix can also be viewed as stochastic
interventions. Instead of replacing A with A∗ , an intervention on features of the network topology replaces A with the members of a class A∗ of n × n adjacency matrices that share the intervention features, stochastically according to some probability distribution gA∗ over A∗ . For example, we might be interested in interventions that constrain the degree distribution of the network, e.g. fixing the maximum degree to be smaller than some number D or ensuring that no
node has degree equal to 0. We might specify gA∗ (A) =

1
|A∗ | I

{A ∈ A∗ }, giving equal weight to

each realization in the class A∗ . Effectively, this kind of intervention sets Vi to v with probability




∗
A∗
A∗
∗
∗
A∗
∗
A∗
P X ∈ x∗ : sT
Y,i (C, x ) = v | W = sX (C) for some A ∈ A , where sX (C) = sX,1 (C), ..., sX,n (C) .
As with the stochastic interventions discussed in the previous section, positivity is a crucial assumption for identifying interventions on A. The support of V ∗ must be the same as the support of
V . If replacing A with A∗ (either deterministically or as a random selection from the class A∗ ) assigns
to unit i a value of V that not observed in the real data for a unit in the same C stratum as i, then
the effect of the intervention that that replaces A with A∗ is not identified for unit i.
We can expand the arguments of sX (·) and sY (·) to include network topology explicitly, in addition
to C, or we can include features of network topology as individual-level covariates. For example, we
could include in Ci the total number of ties for subject i has, a measure of subject i’s centrality in
the network, or a local clustering coefficient for the neighborhood comprised of subject i’s contacts.
Whether or not we can define, identify, and estimate interventions involving these features of network
topology hinges crucially on the positivity assumption.
Both the degree of subject i and local clustering around subject i are local features of network
topology: they depend on T only through subject i and subject i’s immediate contacts. A local
clustering coefficient for node i can be defined as the proportion of potential triangles that include i as
one vertex and that are completed, or the number of pairs of neighbors of i who are connected divided
by the total number of pairs of neighbors of i (Newman, 2009). This measure of triangle completion
captures the extent to which “the friend of my friend is also my friend”: triangle completion is high
whenever two subjects who share a mutual contact are more likely to themselves share a tie than
are two subjects chosen at random from the network. Positivity could hold if, within each level of
C, subjects were observed to have a wide range of degrees and of triangle completion among their

21

contacts.
In contrast with degree and local clustering, network centrality is a node-specific attribute that
nevertheless depends on the entire network topology. It captures the intuitive notion that some nodes
are central and some nodes are fringe in any given network. It can be measured in many different
ways, based, for example, on the number of network paths that intersect node i, on the probability
that a random walk on the network will intersect node i, or on the mean distance between node i and
the other nodes in network (see Chapter 7 of Newman, 2009 for a comprehensive discussion of these
and other centrality measures). Centrality is given by a univariate measure for each node in a network,
but each node’s measure depends crucially on the entire graph. In reality it is not generally possible
to intervene on centrality without altering the entire adjacency matrix A, which is likely to result
in changes to sX (·) and sY (·). Unless it is hypothetically possible to intervene on centrality without
changing any other input into the SEM, in this case sX (·) and sY (·), interventions on centrality cannot
be identified.

4.4

Too many friends, too much influence

The conditions in the theorem on page 16 will be violated for any asymptotic regime in which the degree
of any one or more nodes grows at the same or similar rate to the sample size n. This is problematic
because social networks frequently have a small number of “hubs”–that is, nodes with very high degree
(Newman, 2009), and the occurrence of hubs is a feature of many of the network-generating models
that have been proposed for social networks. When a small number of individuals wield influence over
a significant portion of the rest of the population, two problems arise for statistical inference about
the people connected in the underlying network. First, the number of hubs may grow slowly with n
or it could be fixed and not grow at at all. If the hubs are systematically different from the rest of the
population (e.g. they have different covariates values or the exposure affects them differentially), then
a fixed or slowly growing number of hubs would not allow for consistent inference about this distinct
subpopulation. Second, and more importantly, the sweeping influence of hubs creates dependence
among all of the influenced nodes that undermines inference. Our methods rely on the independence
of Yi and Yj whenever nodes i and j do not share a tie or a mutual alter. When hubs are present, a
significant proportion of nodes will share a connection to one of these hubs, undermining our methods.
We can recover valid inference using our methods if we condition on the hubs, treating them as
features of the background network environment rather than as observations. This results in different
causal effects or statistical estimands, as all of our inference is conditional on the identity and characteristics of the hubs. Imagine a social network comprised of the residents of a city in which a cultural
or political leader is connected to almost all of the other nodes. It may be impossible to disentangle the
influence of this leader, which affects every other node, from other processes simultaneously occurring
among the other residents of the city. It will certainly be impossible to statistically learn about the
hub, as the sample size for the hub subgroup is 1. But it may make sense to consider the hub as a
feature of the city rather than a member of the network. We could then learn about other processes
occurring among the other residents of the city, conditional on the behavior and characteristics of the
leader. For example, we could evaluate the effect of a public health initiative encouraging residents
to talk to their friends about the importance of exercise, but we could not evaluate a similar program

22

targeting the leader’s communication about exercise. It is easy to see that, in the latter case, positivity
would be violated as there would be no control group of subjects unexposed to the leader’s communications. Equally important is the fact that conditioning on the leader is necessary to engender the
independence among units required for our theoretical results to hold.
Practically speaking, for real and finite datasets, this implies that the methods we have proposed are
inappropriate for networks in which the degree is large, compared to n, for one or more nodes. If many
nodes are connected to a significant fraction of other nodes, this problem is intractable. However, if
only a small number of nodes are highly connected we can condition on them to recover approximately
valid inference using our methods for conditional estimands. There is a theoretical tradeoff between the
rate of convergence of our estimators and the order of K relative to n that, in finite samples, becomes
a practical tradeoff between generality and variance. Increasing the number of nodes classified as hubs
(and conditioned on in subsequent analysis) will increase the rate of convergence by decreasing the order
of K for the remaining, non-hub nodes – assuming that the number of hubs remains small compared
to n so that the sample size does not decrease significantly when we exclude hubs from the analysis.
On the other hand, classifying more nodes as hubs results in analyses that are increasingly specific:
conditioning on a single hub may preserve generalizability to other networks (similar cities with similar
leaders), but conditioning on many hubs is likely to limit the generalizability of the resulting inference.

5

Simulations

We conducted a simulation study that evaluated the finite sample and asymptotic behavior of the
TMLE procedure described in Section 3.3. We generated social networks according to the preferential
attachment model (Barabási and Albert, 1999), where the node degree (number of friends) distribution
followed a power law with α = 0.5, and according to the small world network model (Watts and
Strogatz, 1998) with a rewiring probability of 0.1. We generated data with two different types of
dependence: first with dependence due to direct transmission only, and second with both latent variable
dependence and dependence due to direct transmission.
Our simulations mimicked a hypothetical study designed to increase the level of physical activity
in a population comprised of members of a social network. For each community member indexed
by i = 1, . . . , n, the study collected data on i’s baseline covariates, denoted Ci , which included the
indicator of being physically active, denoted P Ai and the network of friends on each subject, Fi .
The exposure or treatment, Xi , was assigned randomly to 25% of the community. For example,
one can imagine a study where treated individuals received various economic incentives to attend
a local gym. The outcome Yi was a binary indicator of maintaining gym membership for a predetermined follow-up period. We assumed that it was of interest to examine and estimate the average
 
of the mean counterfactual outcomes E Ȳ ∗ under various hypothetical interventions g ∗ on such
a community. First, we considered a stochastic intervention g1∗ which assigned each individual to
treatment with a constant probability of 0.35; this differs from the observed allocation of treatment
to 25% of the community members. We also considered a scenario in which the aforementioned
economic incentive was resource constrained and could only be allocated to up to 10% of community
members and estimated the effects of various targeted approaches to allocating the exposure. For

23

example, we considered an intervention g2∗ that targeted only the top 10% most connected members
of the community, as such a targeted intervention would be expected to have a higher impact on the
overall average probability of maintaining gym membership among the community, when compared to
purely random assignment of exposure to 10% of the community. Another hypothetical intervention g3∗
assigned an additional physically active friend to individuals with fewer than 10 friends. Notably, g3∗
can be thought of as intervening on the structure of the social network itself. Finally, we estimated the
combined effect of simultaneously implementing intervention g2∗ and the network-based intervention
g3∗ on the same community; this is g4∗ . For simplicity, this simulation study only reports the expected
outcome under each of these interventions; causal effects defined as contrasts of these interventions
can be easily estimated based on the same methods.
All simulation and estimation was carried out in R language (R Core Team, 2015) with packages
simcausal (Sofrygin et al., 2015) and tmlenet (Sofrygin and van der Laan, 2015). The full R code
for this simulation study is available in a separate github repository1 . The simulations were repeated
for community sizes of n = 500, n = 1, 000 and n = 10, 000. The estimation was repeated by sampling
1, 000 such datasets, conditional on the same network (sampled only once for each sample size). For the
simulations with dependence due to direct transmission, the baseline covariates were independently
and identically distributed. The probability of success for each Yi was a logit-linear function of i’s
exposure Xi (indicator of receiving the economic incentive), the baseline covariates Ci and the three
summary measures of i’s friends exposures and baseline covariates. In particular, we also assumed
that the probability of maintaining gym membership increased on a logit-linear scale as a function of
P
the following network summaries: the total number of i’s friends who were exposed ( j:Aij =1 Xj ), the
P
total number of i’s friends who were physically active at baseline ( j:Aij =1 P Aj ) and the product of
P
P
the two summaries ( j:Aij =1 X j × j:Aij =1 P Aj ). The economic incentive to attend local gym had a
small direct effect on each individual who was not physically active at baseline and no direct effect on
those who were already physically active. However, physically active individuals were more likely to
maintain gym membership over the follow-up period if they had at least one physically active friend
at baseline. We repeated these simulations with the addition of latent variable dependence, which
we introduced by generating unobserved latent variables for each node which affected the node’s own
outcome as well as the outcomes of its friends.
We estimated the expected outcome under interventions g1∗ through g4∗ using the aforementioned
TMLE approach and evaluated its finite sample bias. For the simulations under dependence due
 
to direct transmission, we estimated the marginal parameter E Ȳ ∗ and compared three different
estimators of the asymptotic variance and the coverage of the corresponding confidence intervals.
First, we looked at the naive plug-in i.i.d. estimator (“IID Var ”) for the variance of the influence
curve which treated observations as if they were i.i.d. Second, we used the plug-in variance estimator
based on the efficient influence curve which adjusted for the correlated observations (“dependent IC
Var ”), as previously described in Sofrygin and van der Laan (2015). Finally, we used the parametric
bootstrap variance estimator (“bootstrap Var ”), where each bootstrap iteration sampled n observations
from the previous fit of the likelihood (the estimated exposure model ĝ and the mean outcome model
m̂) and used this bootstrap sample to re-fit the least favorable parametric submodel update . The
1 github.com/osofr/socialnets.sim.study

24

simulation results showing the mean length and coverage of these three CI types are shown in Figure 2
(preferential attachment network) and Figure 3 (small world network). For the simulations with latent
 
variable dependence, we estimated the conditional parameter E Ȳ ∗ and we compared two plug in
variance estimators based on the conditional influence function DC : one that assumes conditionally
i.i.d outcomes (conditional on X and C), which would be true if all dependence were due to direct
transmission but is violated in the presence of latent variable dependence (“IID Var ”), and one that
does not make this assumption (“dependent IC Var ”).
We examined the empirical distribution of the transformed TMLEs, comparing their histogram
estimates against the predicted normal limiting distribution, with the results shown in Figure 4, where
the histogram plots are displayed by sample size (horizontal axis) and the intervention type (vertical
axis). The estimates were first centered at the corresponding true parameter values and then rescaled by their corresponding true standard deviation (SD). We note that our results indicate that the
estimators converge to their normal theoretical limiting distribution, even in networks with power law
node degree distribution, such as the preferential attachment network model, as well as in the densely
connected networks obtained under the small world network model. The results shown in Figure 4
were generated from simulations with dependence due to direct transmission; simulations with latent
variable dependence (not shown) evinced similar approximate normality.
One of the lessons of our simulation study is that by leveraging the structure of the network it might
be possible to achieve a larger overall intervention effect on a population level (Harling et al., 2016).
For example, the results in the left panels of Figures 2 and 3 show that by targeting the exposure
assignment to highly connected and physically active individuals, intervention g2∗ increases the mean
probability of sustaining gym membership compared to the similar level of un-targeted coverage of the
exposure. We also demonstrated the feasibility of estimating effects of interventions on the observed
network structure itself, such as intervention g3∗ , which can be also combined with economic incentives,
as it was mimicked by our hypothetical intervention g2∗ + g3∗ . These combined interventions could be
particularly useful in resource constrained environments, since they may result in larger community
level effects at the lower coverage of the exposure assignment.
We evaluated the performance of the different approaches to variance estimation described above,
as measured by the coverage of the 95% CIs. Our results from simulations with dependence due to
direct transmission show that conducting inference while ignoring the nature of the dependence in such
datasets generally results in anticonservative variance estimates and under-coverage of CIs, which can
be as low as 50% even for very large sample sizes (“IID Var ” in the right panels of Figures 2 and 3).
The CIs based on the dependent variance estimates (“dependent IC Var ” in the right panel of the same
figures) obtain nearly nominal coverage of 95% for large enough sample sizes, but can suffer in smaller
sample sizes due to lack of asymptotic normality and near-positivity violations. Notably, the CIs based
on the parametric bootstrap variance estimates provide the most robust coverage for smaller sample
sizes, while attaining the nominal 95% coverage in large sample sizes for nearly all of the simulation
scenarios (“bootstrap Var ” in the right panel of the same figures). The apparent robustness of the
parametric bootstrap method for inference in small sample sizes, even as low as n = 500, was one of
the surprising finding of this simulation study. Future work will explore the assumptions under which
this parametric bootstrap works and its sensitivity towards violations of those assumptions. Similarly,

25

CI.type

dependent IC Var

bootstrap Var

iid Var

g1∗ (random 35%)

g3∗ (network intervention)

N: 500

N: 500

g2∗ (dynamic intervention)

g2∗ + g3∗
g1∗ (random 35%)

g3∗ (network intervention)

N: 1000

N: 1000

g2∗ (dynamic intervention)

g2∗ + g3∗
g1∗ (random 35%)

g3∗ (network intervention)

N: 10000

N: 10000

g2∗ (dynamic intervention)

g2∗ + g3∗
0.00

0.05

0.10

0.15

Mean estimate & 95% CI length

0.20

0.7

0.8

0.9

Coverage

Figure 2: Mean 95% CI length (left panel) and coverage (right panel) for the TMLE in preferential
attachment network with dependence due to direct transmission, by sample size, intervention and CI
type. Results are shown for the estimates of the average expected outcome under four hypothetical
interventions (g1∗ , g2∗ , g3∗ and g2∗ + g3∗ ).
in the simulations with latent variable dependence the variance estimates that assume conditionally
i.i.d. outcomes, i.e. that dependence may be due to direct transmission but not to latent variables,
are anti-conservative.

6

Conclusion

We proposed new methods that allow for causal and statistical inference using observations sampled
from members of a single interconnected social network when the observations evince dependence
due to network ties. In contrast to existing methods, our methods do not require randomization of an
exogenous treatment and they have proven performance under asymptotic regimes in which the number
of network ties grows (slowly) with sample size. In future work we plan to address a key limitation
of the present proposal, namely the assumption that the network is observed fully and without error.
We also plan to develop data-adaptive methods for estimating the summary measures sX and sY , as
it may be unreasonable to expect these to be known a priori. Finally, we plan to develop estimating
algorithms for longitudinal settings; the influence function and asymptotic results for these settings
are straightforward extensions of the results presented here, but estimation can be challenging.

26

CI.type

dependent IC Var

bootstrap Var

iid Var

g1∗ (random 35%)

g3∗ (network intervention)

N: 500

N: 500

g2∗ (dynamic intervention)

g2∗ + g3∗
g1∗ (random 35%)

g3∗ (network intervention)

N: 1000

N: 1000

g2∗ (dynamic intervention)

g2∗ + g3∗
g1∗ (random 35%)

g3∗ (network intervention)

N: 10000

N: 10000

g2∗ (dynamic intervention)

g2∗ + g3∗
0.0

0.2

0.4

0.6

Mean estimate & 95% CI length

0.8

0.6

0.7

0.8

0.9

Coverage

Figure 3: Mean 95% CI length (left panel) and coverage (right panel) for the TMLE in small world
network with dependence due to direct transmission, by sample size, intervention and CI type. Results
are shown for the estimates of the average expected outcome under four hypothetical interventions (g1∗ ,
g2∗ , g3∗ and g2∗ + g3∗ ).

Figure 4: Comparing re-scaled empirical TMLE distributions (black) to their theoretical normal limit
(red) with varying sample size (x-axis) and intervention type (y-axis). TMLEs were centered at the
truth and then re-scaled by true SD. Results shown for the preferential attachment network (left) and
the small world network (right).

27

CI.type

dependent IC Var

iid Var

g1∗ (random 35%)

g3∗ (network intervention)

N: 500

N: 500

g2∗ (dynamic intervention)

g2∗ + g3∗
g1∗ (random 35%)

g3∗ (network intervention)

N: 1000

N: 1000

g2∗ (dynamic intervention)

g2∗ + g3∗
g1∗ (random 35%)

g3∗ (network intervention)

N: 10000

N: 10000

g2∗ (dynamic intervention)

g2∗ + g3∗
0.00

0.05

0.10

0.15

Mean estimate & 95% CI length

0.20

0.75

0.80

0.85

0.90

0.95

Coverage

Figure 5: Mean 95% CI length (left panel) and coverage (right panel) for the TMLE in preferential
attachment network with latent variable dependence, by sample size, intervention and CI type. Results
are shown for the estimates of the average expected outcome under four hypothetical interventions (g1∗ ,
g2∗ , g3∗ and g2∗ + g3∗ ).

28

CI.type

dependent IC Var

iid Var

g1∗ (random 35%)

g3∗ (network intervention)

N: 500

N: 500

g2∗ (dynamic intervention)

g2∗ + g3∗
g1∗ (random 35%)

g3∗ (network intervention)

N: 1000

N: 1000

g2∗ (dynamic intervention)

g2∗ + g3∗
g1∗ (random 35%)

g3∗ (network intervention)

N: 10000

N: 10000

g2∗ (dynamic intervention)

g2∗ + g3∗
0.0

0.2

0.4

0.6

Mean estimate & 95% CI length

0.6

0.7

0.8

0.9

Coverage

Figure 6: Mean 95% CI length (left panel) and coverage (right panel) for the TMLE in small world
network with latent variable dependence, by sample size, intervention and CI type. Results are shown
for the estimates of the average expected outcome under four hypothetical interventions (g1∗ , g2∗ , g3∗
and g2∗ + g3∗ ).

29

Acknowledgements
Elizabeth Ogburn was supported by ONR grant N000141512343. Oleg Sofrygin and Mark van der
Laan were supported by NIH grant R01 AI074345-07.

Appendix: Proof of Theorem 1
Regularity conditions
For a real-valued function c 7→ f (c), let the L2 (P )-norm of f (c) be denoted by kf k = E[f (C)2 ]1/2 .
Define Mm and Mh̃ as the classes of possible functions that can be used for estimating the two nuisance
parameters m and h̃ ≡ h̄x∗ /h̄, respectively. Note that a model for g plus the empirical distribution of
covariates C determines h̃. Equivalent assumptions could be stated in terms of g instead of h̃, but we
focus on h̃ because that is the functional of g and C that we model in our estimating procedure. Assume
ˆ ∗ /h̄
ˆ ∈ M with probability
that the TMLE update m̂ ∈ M with probability 1 and assume that h̄
ˆ

m

x

h̃

1. Finally, define the following dissimilarity measure on the cartesian product of F ≡ Mm × Mh̃ :





d (h, m) , h̃, m̃ = max sup | h − h̃ | (v), sup | m − m̃ | (v) .
v∈V

v∈V

The following are the regularity conditions required for Theorem 1, i.e. for asymptotic normality
of the TMLE ψ̂ ∗ .
Uniform consistency: Assume that
d





ˆ ∗ /h̄
ˆ , m̂ , h̄ ∗ /h̄, m → 0
h̄
x
ˆ
x

in probability as n → ∞. Note that this assumption is only needed for proving the asymptotic
equicontinuity of our process; it is not needed for proofs of relevant convergence rates for the
second order terms.
Bounded entropy integral: Assume that there exists some η > 0, so that

´ηp
log (N (, F, d))d <
0

∞, where N (, F, d) is the number of balls of size  w.r.t. metric d needed to cover F.
Universal bound: Assume supf ∈F ,O | f | (O) < ∞, where the supremum of O is over a set that
contains O with probability one. This assumption will typically be a consequence of the choosing
a specific function class F that satisfies the above entropy condition.
Positivity: Assume
sup
v∈V

h̄x∗ (v)
< ∞.
h̄(v)

ˆ − h̄ =
Consistency and rates for estimators of nuisance parameters: Assume that km̂ − mk h̄


−1/2
oP (Cn )
. Note that this rate is achievable if, for example, estimation of h̄ relies on some

30

−1/4

pre-specified parametric model, or if both h̄ and m are estimated at rate Cn

.

Rate of the second order term: Assume that
)
!
ˆ ( ˆ
 p 
h̄x∗
h̄x∗
−
(m̂ˆ − m)(v)h̄(v)dµ(v) = oP 1/ Cn .
Rn1 ≡ −
ˆ
h̄
v
h̄
Note that this condition is provided here purely for the sake of completeness, since it will satisfied
ˆ − h̄ . This follows from
based on the previously assumed rates of convergence for km̂ − mk h̄
the fact that the parametric TMLE update step m̂ˆ of m̂ will have a negligible effect on the rate
of convergence of the initial estimator m̂, that is, m̂ˆ will converge at “nearly” the same rate as
m̂.
Limited connectivity and limited dependence of Y,X and C: Let Kmax,n = maxi {Ki } for a
2
fixed network with n nodes. Assume that Kmax,n
/n converges to 0 in probability as n → ∞.

Bounded fourth moments: E[fi (O)4 ] < ∞, where fi (O) is the contribution of the ith observation
to the estimator and is defined below.
A key condition is consistency and rates for estimators of nuisance parameters. This condition will
be satisfied, for example, if both nuisance models are parametric and if one of the two is correctly
1/4

specified, or if both models converge to the truth at rate Cn . It can in fact be weakened, but for
a more general discussion and the corresponding technical conditions we refer to the Appendix of
van der Laan (2014). With the exception of the rates of convergence, the more general conditions for
asymptotic normality of the TMLE presented in that paper apply to our setting as well.

7

Overview of the proof of Theorem 1

√
We want to show that
Cn (ψ̂ − ψ) converges in law to a Normal limit as n goes to infinity for some
q
√
√
√
√
2
rate Cn such that n/ (Kmax (n)) ≤ C n ≤ n, where the rate Cn is the order of the variance
of the sum of the first-order linear approximation of (ψ̂ − ψ).
Broadly, the proof has two parts: First, we require that the second order terms in the expansion
√
of ψ̂ − ψ are stochastically less than 1/ Cn , that is that
ψˆn − ψ

=

n
 p 
1X
{fi (O) − E[fi (O)]} + op 1/ Cn ,
n i=1

where fi (O) is the contribution of the ith observation to the estimator. Specifically, for our influence
function
D(o)

=

n
n
n
X
1X
1 X q̄x∗ (vi )
E [m (Vi∗ ) | Cj = cj ] − ψ +
{yi − m (vi )} ,
n i=1
n i=1 q̄(vi )
j=1

31

the contribution of the ith observation is
fi (o)

=

n
X

E [m (Vi∗ ) | Cj = cj ] +

i=j

q̄x∗ (vi )
{yi − m (vi )} .
q̄(vi )

Then proving asymptotic normality of the TMLE amounts to the asymptotic analysis of the sum
Pn
i=1 {fi (O) − E[fi (O)]}, and the second part of the proof establishes that the first order terms
√ Pn
√
converge to a normal distribution when scaled by Cn , that is that Cn i=1 {fi (O) − E[fi (O)]} →d
N (0, σ 2 ) for some finite σ 2 .

√
The proof that the second order terms are stochastically less than 1/ Cn is an extension of the

empirical process theory of Van Der Vaart and Wellner (1996) and follows the same format as the
proof in van der Laan (2014). Indeed, the proof offered by van der Laan (2014) holds immediately after
√
√
replacing the rate or scaling factor n with Cn throughout. Only one step in the van der Laan (2014)
proof relies on the network topology, which is the major difference between the setting in that paper,
where the number of network connections is fixed and bounded as n goes to infinity, and the present
setting: the proof requires bounding the Orlicz norms of several empirical processes corresponding to
p

components of the influence function for ψ, and a key step is bounding the expectation of E [|Xn (f )| ]
, where Xn (f ) is the stochastic process that describes the difference between the empirical (indexed by
n) and the true distribution functions of a component of the influence function for ψ. This step relies
on a combinatorial argument about nature of overlapping friend groups in the underlying network, and
the argument for the case of growing Ki is subsumed by the argument for fixed K in van der Laan
(2014).
The proof that the first order terms converge to a normal distribution requires a central limit
theorem for dependent data with growing and possibly irregularly sized dependency neighborhoods,
where a dependency neighborhood for unit i is a collection of observations on which the observations
for unit i may be dependent. We prove such a CLT in Lemmas 1 and 2. In next section we use the
CLT for growing and irregular dependency neighborhoods, along with an orthogonal decomposition of
the first order terms, to prove the remainder of Theorem 1.

Central limit theorem for first order terms
Pn
Proving asymptotic normality of the TMLE amounts to the asymptotic analysis of the sum i=1 {fi (O) − E[fi (O)]}.
Pn
As a start, decompose i=1 {fi (O) − E[fi (O)]} into a sum of three orthogonal components:
fY,i (Y, X, C)

=

fi (O) − E [fi (O) | X, C] ,

fX,i (X, C)

=

E[fi (O) | X, C] − E[fi (O) | C], and

fC,i (C)

=

E[fi (O) | C] − E[fi (O)].

Note that
fi (O) − E[fi (O)] = fY,i (Y, X, C) + fX,i (X, C) + fC,i (C)

32

and with slight abuse of notation we will also write fY,i (O), fX,i (O) and fC,i (O). Let fY (O) =
Pn
Pn
Pn
i=1 fY,i (O), fX (O) =
i=1 fX,i (O) and fC (O) =
i=1 fC,i (O). For i = 1, . . . , n, let

Z Y,i

=

Z X,i

=

Z C,i

=

fY,i (Y, X, C)
Pn
V ar( i=1 fY,i (Y, X, C))
fX,i (X, C)
p
Pn
V ar( i=1 fX,i (X, C))
fC,i (C)
p
.
Pn
V ar( i=1 fC,i (C))
p

and

0
ZY,i

=

0
ZX,i

=

fY,i (Y, X, C) |(X, C)
Pn
V ar( i=1 fY,i (Y, X, C) |(X, C) )
fX,i (X, C) |C
p
Pn
V ar( i=1 fX,i (X, C) |C )
p

0
conditions fY,i (O) on (X, C) and
We use the prime to denote conditional random variables: ZY,i
0
conditions fX,i (O) on C and
rescales it by the standard error of fY (O)| (X, C). Similarly, ZX,i

rescales it by the standard error of fX (O)|C. Let
2
σnY

(x, c) = V ar

n
X

!
fY,i (Y, x, c) |(X = x, C = c)

i=1

 2

2
σnY
= EPX,C σnY
(X, C) ,

2
σnX
(c)

= V ar

n
X

!
fX,i (X, c) |C = c

i=1

 2

2
σnX
= EPC σnX
(C) ,
and
2
σnC

= V ar

n
X

!
fC,i (C) .

i=1

Pn
Pn
2
2
Note that by the law of total variance σnX
= V ar( i=1 fX,i (X, C)) and σnY
= V ar ( i=1 fY,i (Y, X, C)).
Pn
P
P
P
n
n
n
0
0
0
0
, ZnX
denote i=1 ZX,i
, ZnY denote i=1 ZY,i , ZnX denote i=1 ZX,i ,
Let ZnY
denote i=1 ZY,i
Pn
and ZnC denote i=1 ZC,i . We will establish convergence in distribution of each of the three terms
0

0

separately. Because ZnY and ZnX converge to distributions that do not depend on their conditioning
events, conditional convergence in distribution implies convergence of ZnY and ZnX to the same
limiting distributions. Since fY (O),fX (O), and fC (O) are orthogonal by construction, the variance
of the limiting distribution of their sum is the sum of their marginal variances. If the three processes
converge at the same rate the limiting variance will be the sum of the variances of the three processes.
33

However, the three terms may converge at different rates, in which case the limiting distribution of
ψ̂ − ψ will be given by the limiting distribution of the term(s) with the slowest rate of convergence.
0
0
In order to show that ZnX
, ZnY
, and ZnC all converge in distribution to a N (0, 1) random variable,

we can use three separate applications of the central limit theorem given in Lemma 1, which is based
on Stein’s method.
Stein’s method (Stein, 1972) quantifies the error in approximating a sample average with a normal
distribution. (For an introduction to Stein’s method see Ross, 2011.) Stein’s method has been used to
prove CLTs for dependent data with dependence structure given by dependency neighborhoods (Chen
and Shao, 2004): the dependency neighborhood for observation i is a set of indices Di such that
observation i is independent of observation j, for any j ∈
/ Di . Conditionally on C, fX,i and fX,j are
independent for any nodes i and j such that Aij = 0 and there is no k with Aik = Ajk = 1, that is
for any nodes that do not share a tie or have any mutual network contacts. The same is true for fY,i
and fY,j conditional on X and C and for fC,i and fC,j . Thus the three collections of random variables
0
0
0
0
ZX,1
, ..., ZX,n
, ZY,1
, ..., ZY,n
, and ZC,1 , ..., ZC,n each has a dependency neighborhood structure with

Di = i ∪ {j : Aij = 1} ∪ {k : Ajk = 1 for j : Aij = 1}, that is the “friends” and “friends of friends”
of node i. Define the indicators R(i, j) for any (i, j) ∈ {1, . . . , n}2 to be an indicator of dependence
between ZX,i and ZX,j , R(i, j) = 1 iff j ∈ Di or, equivalently, if i ∈ Dj . For any i ∈ {1, . . . , n} the set
0
0
{ZX,j
: (R(i, j) = 1, j ∈ {1, . . . , n})} forms the dependency neighborhood of ZX,i
and the collection
0
0
{ZX,j
: (R(i, j) = 0, j ∈ {1, . . . , n})} is independent of ZX,i
. The same logic applies to defining the
0
0
dependency neighborhoods for ZY,1
, ..., ZY,n
conditional on X and C, and for ZC,1 , ..., ZC,n based on

(unconditional) independence of each fC,i (O) and fC,j (O), as determined by the network structure
and the distributional assumptions made for the baseline covariates C.
0
Applied to ZnX
, Stein’s method provides the following upper bound

0
d(ZnX
, Z) ≤

n
X
X

0
0
0
E ZX,i
ZX,j
ZX,k

i=1 j,k∈Di

v


r u
n X
u
X
2u
0 Z 0 ,
tV ar 
+
ZX,i
X,j
π
i=1
j∈Di

where Z ∼ N (0, 1) and d(·, ·) is the Wasserstein distance metric (Vallender, 1974).
0
In order to show that ZnX
converges in distribution to Z, we must show that the righthand side

of the inequality converges to zero as n goes to infinity. We will first show that this convergence holds
when Ki = |Fi | = Kmax (n) for all i, that is when all nodes have the same number of ties. We will
then show that removing any tie from the network preserves an upper bound on the righthand side of
the inequality. This completes our proof that for any network such that Ki ≤ Kmax (n) for all i and
2
Kmax
(n)
0
converges to zero as n goes to infinity, ZnX
converges in distribution to a standard normal disn
tribution. The same argument applied to Z nC proves that it has a Normal limiting distributions as well.

1. [Applying Stein’s Method to the dependent sum] Consider a network of nodes given by adjacency
matrix A. Let U1 , ..., Un be bounded mean-zero random variables with finite fourth moments and with
34

dependency neighborhoods Di = i ∪ {j : Aij = 1} ∪ {k : Ajk = 1 for j : Aij
= 1}, and let Ki be the
P
D
Ui
2
degree of node i. If Ki = Kmax (n) for all i and Kmax (n) /n → 0, then √
→ N (0, 1).
P
var(

Proof. [Proof of Lemma 1] Let Ui0 = √

Ui
.
P
var( Ui )

Ui )

Application of Stein’s method often involves defining

0

the so-called “Stein coupling” (W, W , G) (Fang, 2011; Fang et al., 2015). Consider the following sum
Pn
of dependent variables W = i=1 Ui0 . Define a discrete random variable I distributed uniformly over
Pn
{1, . . . , n} and define another random variable W 0 = (W − j=1 R(I, j)Uj0 ). Finally, define G = −nUI0
and note that (W, W 0 , G) forms a Stein coupling \citep{fang2011thesis,fang2015rates}. We also let
PN
D = (W 0 − W ) = − j=1 R(I, j)U 0 j . This Stein coupling allows us then to derive the upper bound

d(W, Z) ≤

n
X

X

E U 0 i Uj0 Uk0

i=1 j,k∈Di

v


r u
n X
u
X
2u
tV ar 
+
Ui0 Uj0 ,
π
i=1

(9)

j∈Di

as shown in Ross (2011). We will now show that, for any network structure,
n
X

v


r u
n X
u
X
X
2
u
tV ar 
E |U 0 i U 0 j Uk0 | +
Ui0 Uj0 
π
i=1

i=1 j,k∈Di

j∈Di


=



P
 i,j,k R(i, j)R(i, k) 
O h
i3/2  .
P
i,j R(i, j)

The righthand side of the above equation is equal to

q

(10)

(Kmax (n))2
n

under the assumption of Kmax (n)

ties for each node i = {1, . . . , n}. By assumption, we also have that

Kmax (n)
√
n

converges to zero as n

goes to infinity, and therefore if we can show equation (10) we have proved that √

D
Ui
→
P
var( Ui )
P

N (0, 1).

Consider the term
n
X

X

i=1 j,k∈Di

E |U 0 i U 0 j Uk0 | =

1
P
var( Ui )3/2

2 


X
Uk 
.
E Ui 




i=1
j∈Di

n
X







h
i3/2 
P
P
P
By the assumption of bounded 4th moments, var( Ui )3/2 = O
R(i,
j)
, that is, var( Ui )
i,j
P
stabilizes to a constant when scaled by i,j R(i, j). Using the fact that each |Ui | is bounded we get
2 


X
E Ui 
Uj 




i=1
j∈Di


n X

X
≤ M
R(i, j)R(i, k)


i=1
j,k
X
= M
R(i, j)R(i, k),
N
X







i,j,k

35

for some positive constant M < ∞. Combining the above expressions, we get
n
X




P

 i,j,k R(i, j)R(i, k) 
E |U 0 i U 0 j Uk0 | = O  h
i3/2  .
P
i=1 j,k∈Di
R(i,
j)
i,j
X

Now consider the second term:
r
P

v


u
n P
V
ar
U
U
n
u
i
j
i=1
j∈Di
XX
u
P
tV ar 
Ui0 Uj0  =
.
var( Ui )2
i=1
j∈Di

P

i,j

R(i, j) terms in

Pn

P

U i U j , and the number of terms Uk Ul with which Ui Uj has
P

P
P
P
n P
non-zero covariance is |Di ∪ Dj | ≤ k R(i, k)+ k R(i, k), so V ar
i=1
j∈Di U i U j ≤ M
i,j R(i, j)
k R(i, k)
P

P

P
P
n
for some finite M . Therefore V ar
= O
Ui )2 =
i=1
j∈Di U i U j
i,j,k R(i, j)R(i, k) . V ar(
h

i2
P
O
, so the second term is of smaller order than the first term. Therefore we have
i,j R(i, j)

There are

i=1

j∈Di

P

only to consider the first term and we have completed the proof.
2. [Bound goes to zero when Ki ≤ Kmax (n) for all i] Convergence to zero of the righthand side of
Equation (9) is preserved under the removal of ties and holds as long as Ki ≤ Kmax (n) for all i and
2
Kmax
(n)
n

converges to zero as n goes to infinity.

Proof. [Proof of Lemma 2] Consider a sequence of networks with n going to infinity such that the
righthand side of Equation (9) converges to 0, i.e.
n
X
X

E Ui0 Uj0 Uk0

i=1 j,k∈Di

v


r u
n X
u
X
2u
tV ar 
Ui0 Uj0 
+
π
i=1

→ 0.

j∈Di

Because the second term is of the same or smaller order than the first, we only have to consider
Pn P
the first term. For this sequence of networks, define An = i=1 j,k∈Di E Ui0 Uj0 Uk0 . Removing a
single tie from the underlying network has the effect of rendering independent some pairs that were
previously dependent; We now consider the effect of rendering a single dependent pair independent
but otherwise leaving the distributions of the random variables the same. Suppose the pair rendered
independent is (l, m). Define a new sequence of networks with n going to infinity to be identical to
the previous sequence but with pair (l, m) independent, and let A0n be the first term in the righthand
side of Equation (9) for this new sequence. Then
A0n = An − 2

X

0
E |Ul0 Um
Uk0 |

k∈Dl ∪Dm

which is bounded above by An .
0
0
This completes the proof that ZnX
, ZnY
, and Z nC have Normal limiting distributions.

36

0
3. [Conditional CLT implies marginal CLT]ZnX
converges to Normal distribution after marginalizing
0
over C (but conditioning on the network as captured by the adjacency matrix A) and ZnY
converges to

Normal distribution after marginalizing over (X, C). That is, ZnX and ZnY both converge to Normal
distributions.
0
Proof. [Proof of Lemma 3] For illustration consider ZnX
=

Pn

i=1

0
Z2,i
, where

q
0
2 (C)
ZX,i
= (fX,i (X, C) |C ) / σnX
and note that the proof of the convergence of ZnY is nearly identical. The conditional CLT results
from Lemma 1 show that
"
P

0
[ZnX

≤ x |C = c ] = P

!
#
N
X
fX,i (X, c)
p
≤ x |C = c
2 (c)
σnX
i=1

converges to Φ(x) for each x and almost every c, where Φ is the cumulative distribution function of
the standard Normal random variable and C is a given sequence (Ci : i = 1, . . . , n). Let PC denote
the distribution of C. Then
"
P (ZnX ≤ x) ≡ P
ˆ

N
X
fX,i (X, C)
p
≤x
2
σnX
i=1

!#

0
P (ZnX
≤ x|C = c)dPC (c).

=
c

0
For a given x, the dominated convergence theorem is now applied with fn (c) = P (ZnX
≤ x|C = c)

and the limit given by f (c) = Φ(x) = m, where m is some constant that doesn’t depend on c. From
the previous conditional CLT result it follows that fn (c) converges to f (c) pointwise for each c. The
´
next step is to find an integrable function g, such that fn < g and g(c)dPC (c) < ∞. The proof is
then completed by choosing g = 1.
We have now shown that ZnY , ZnX , and ZnC are asymptotically Normally distributed. We now
show that the sum of the three processes converges in distribution to a Normal random variable.
Consider three cases: (1) the three processes have the same rate of marginal convergence in distribution,
(2) one of the three processes converges faster than the other two, and (3) two of the processes converge
faster than the third. In all three cases the rate of convergence for the sum will be the slowest of the
three marginal rates. In case (3), the limiting distribution of the sum is determined entirely by the one
process that converges with a slower rate than the other two: the other two processes will converge to
constants (specifically to their expected values of 0) when standardized by the slower rate; Slutsky’s
theorem concludes the proof. We focus on case (1) below; case (2) follows immediately by applying
the proof below to the two processes that converge at the same slower rate and applying Slutsky’s to
the third, faster converging process.
For convenience, in order to show that the sum of the three dependent processes also converges to
Normal, define
2
2
2
Cn∗ := σnY
+ σnX
+ σnC
.

37

Note that Cn∗ is related to Cn as follows: Cn = O(n2 /Cn∗ ).
4. [CLT for the sum of the three orthogonal processes]If all three processes have the same marginal
rate of convergence, then
1
p
(fY (Y, X, C) + fX (X,C) + fC (C)) → N (0, 1).
Cn∗
Proof. [Proof of Lemma 4] Without the loss of generality, we prove that ZnX + ZnC → N (0, 2) and
note that the general result for (ZnY + ZnX + ZnC ) follows by applying a similar set of arguments.
Consider the following random vector (ZnX , ZnC ) taking values in IR2 . Let Fn (x1 , x2 ) ≡ P (ZnX ≤
x1 , ZnC ≤ x2 ), where (x1 , x2 ) ∈ IR2 . Let Φ2 (x1 , x2 ) ≡ P (ZX ≤ x1 )P (ZC ≤ x2 ), for ZX ∼ N (0, 1)
and ZC ∼ N (0, 1), that is, Φ2 (x1 , x2 ) defines the CDF of the bivariate standard normal distribution,
for (x1 , x2 ) ∈ IR2 . The goal is to show that Fn (x1 , x2 ) → Φ2 (x1 , x2 ), for any (x1 , x2 ) ∈ IR2 . The
convergence in distribution for ZnX + ZnC will follow by applying the Cramer and Wold Theorem
(1936).
Note that
P (ZnX ≤ x1 , ZnC ≤ x2 )
=P (ZnX ≤ x1 |ZnC ≤ x2 )P (ZnC ≤ x2 ).
First, from the previous application of Stein’s method, we have that
P (ZnC ≤ x2 ) → Φ(x2 ),
where Φ(x2 ) ≡ P (ZC ≤ x2 ), ZC ∼ N (0, 1) and x2 ∈ IR2 . Also note that
P (ZnX ≤ x1 |ZnC ≤ x2 )
X
=
P (ZnX ≤ x1 |C = c )P (C = c |ZnC ≤ x2 ),
c∈C

where C denotes the support of C, ZnX = √1 ∗ fX (X,C), ZnC = √1 ∗ fC (C) and
Cn

P (C = c |ZnC ≤ x2 ) =

Cn

p 
P (C = c)I( 1/ Cn∗ fC (c) ≤ x2 )
p 
.
P ( 1/ Cn∗ fC (c) ≤ x2 )

By another application of Stein’s method, it was shown that
P (ZnX ≤ x1 |C = c ) → Φ(x2 ),
for any realization of c ∈ C. That is, we’ve shown that the limiting distribution of ZnX conditional
on C = c, does not itself depend on the conditioning event C = c. Applying Lemma 3, we finally
conclude that Fn (x1 , x2 ) → Φ2 (x1 , x2 ), for any (x1 , x2 ) ∈ IR2 and the result follows.

38

Variance estimation
The estimate of the variance of the TMLE ψ̂ can be obtained from the sum, scaled by 1/n2 , of the
three plug-in estimators of
2
σnY

=

X

E(fY,i (O)fY,j (O))

i,j
2
σnX

=

X

E(fX,i (O)fX,j (O))

i,j
2
σnC

=

X

E(fC,i (O)fC,j (O)).

i,j

Alternatively, one can estimate the variance from a single plug-in estimator
1 X
E(fi (O)fj (O)).
n2 i,j
Note that contribution to these variances of any pair i, j not in each others dependency neighborhoods
will be 0. Therefore, it is acceptable to sum only over pairs i, j sharing a tie or a mutual contact in
√
the underlying network. Finally, note that we do not need to know the true rate of convergence Cn
to obtain a valid estimate of the C.I. for ψ; this rate is captured by the number of non-zero terms in
the variance sums.

References
Ali, M. M. and Dwyer, D. S. “Social network effects in alcohol consumption among adolescents.”
Addictive behaviors, 35(4):337–342 (2010).
Aronow, P. M. and Samii, C. “Estimating average causal effects under general interference.” Technical
report, Yale University (2013).
Athey, S., Eckles, D., and Imbens, G. W. “Exact P-values for Network Interference*.” Journal of the
American Statistical Association, (just-accepted) (2016).
Barabási, A.-L. and Albert, R. “Emergence of scaling in random networks.” science, 286(5439):509–512
(1999).
Besag, J. “On spatial-temporal models and Markov fields.” In Transactions of the Seventh Prague
Conference on Information Theory, Statistical Decision Functions, and Random Processes, 47–55.
Springer (1974).
Bowers, J., M, F. M., and C, P. “Reasoning about interference between units: A general framework.”
Political Analysis, 21:97–124 (2013).
Cacioppo, J. T., Fowler, J. H., and Christakis, N. A. “Alone in the crowd: the structure and spread of
loneliness in a large social network.” Journal of personality and social psychology, 97(6):977 (2009).

39

Chen, L. H. and Shao, Q.-M. “Normal approximation under local dependence.” The Annals of Probability, 32(3):1985–2028 (2004).
Christakis, N. and Fowler, J. “The spread of obesity in a large social network over 32 years.” New
England Journal of Medicine, 357(4):370–379 (2007).
—. “The collective dynamics of smoking in a large social network.” New England journal of medicine,
358(21):2249–2258 (2008).
—. “Social network sensors for early detection of contagious outbreaks.” PloS one, 5(9):e12948 (2010).
Clauset, A., Shalizi, C. R., and Newman, M. E. “Power-law distributions in empirical data.” SIAM
review , 51(4):661–703 (2009).
Cohen-Cole, E. and Fletcher, J. “Is obesity contagious? Social networks vs. environmental factors in
the obesity epidemic.” Journal of Health Economics, 27(5):1382–1387 (2008).
Diaconis, P. and Janson, S.

“Graph limits and exchangeable random graphs.”

arXiv preprint

arXiv:0712.2749 (2007).
Eckles, D., Karrer, B., and Ugander, J. “Design and analysis of experiments in networks: Reducing
bias from interference.” arXiv preprint arXiv:1404.7530 (2014).
Fang, X. “Multivariate, combinatorial and discretized normal approximations by Stein’s method.”
Ph.D. thesis (2011).
Fang, X., Röllin, A., et al. “Rates of convergence for multivariate normal approximation with applications to dense graphs and doubly indexed permutation statistics.” Bernoulli , 21(4):2157–2189
(2015).
Goetzke, F. “Network effects in public transit use: evidence from a spatially autoregressive mode
choice model for New York.” Urban Studies, 45(2):407–417 (2008).
Graham, B., Imbens, G., and Ridder, G. “Measuring the Effects of Segregation in the Presence of
Social Spillovers: A Nonparametric Approach.” Technical report, National Bureau of Economic
Research (2010).
Halloran, M. and Hudgens, M. “Causal Inference for Vaccine Effects on Infectiousness.” The University
of North Carolina at Chapel Hill Department of Biostatistics Technical Report Series, 20 (2011).
Halloran, M. and Struchiner, C. “Causal inference in infectious diseases.” Epidemiology, 142–151
(1995).
Haneuse, S. and Rotnitzky, A. “Estimation of the effect of interventions that modify the received
treatment.” Statistics in medicine, 32(30):5260–5277 (2013).
Harling, G., Wang, R., Onnela, J.-P., and DeGruttola, V. “Leveraging Contact Network Structure in
the Design of Cluster Randomized Trials.” Harvard University Biostatistics Working Paper Series,
(Working Paper 199) (2016).
40

Hong, G. and Raudenbush, S. “Evaluating Kindergarten Retention Policy.” Journal of the American
Statistical Association, 101(475):901–910 (2006).
—. “Causal inference for time-varying instructional treatments.” Journal of Educational and Behavioral
Statistics, 33(3):333–362 (2008).
Hudgens, M. and Halloran, M. “Toward causal inference with interference.” Journal of the American
Statistical Association, 103(482):832–842 (2008).
Knoke, D. and Yang, S. Social network analysis, volume 154. Sage (2008).
Lauritzen, S. L. and Richardson, T. S. “Chain graph models and their causal interpretations.” Journal
of the Royal Statistical Society: Series B , 64(3):321–348 (2002).
Lee, L.-F. “Asymptotic Distributions of Quasi-Maximum Likelihood Estimators for Spatial Autoregressive Models.” Econometrica, 72(6):1899–1925 (2004).
Lin, X. “Peer effects and student academic achievement: an application of spatial autoregressive model
with group unobservables.” Unpublished manuscript, Ohio State University (2005).
Lovász, L. Large networks and graph limits, volume 60. American Mathematical Soc. (2012).
Lyons, R. “The spread of evidence-poor medicine via flawed social-network analysis.” Statistics,
Politics, and Policy, 2(1) (2011).
Madan, A., Moturu, S. T., Lazer, D., and Pentland, A. S. “Social sensing: obesity, unhealthy eating
and exercise in face-to-face networks.” In Wireless Health 2010 , 104–110. ACM (2010).
Moreno, J. L. “Sociometry in relation to other social sciences.” Sociometry, 1(1/2):206–219 (1937).
Muñoz, I. D. and van der Laan, M. “Population intervention causal effects based on stochastic interventions.” Biometrics, 68(2):541–549 (2012).
Newman, M. Networks: an introduction. Oxford: Oxford University Press (2009).
Newman, M. E. and Park, J. “Why social networks are different from other types of networks.” Physical
Review E , 68(3):036122 (2003).
Ogburn, E. and VanderWeele, T. J. “Causal diagrams for interference.” Technical report, Harvard
University (2013).
Ogburn, E. L. and VanderWeele, T. J. “Vaccines, Contagion, and Social Networks.” arXiv preprint
arXiv:1403.1241 (2014).
O’Malley, J. A. and Marsden, P. V. “The analysis of social networks.” Health services and outcomes
research methodology, 8(4):222–269 (2008).
Pearl, J. “Causal diagrams for empirical research.” Biometrika, 82(4):669–688 (1995).
—. Causality: models, reasoning and inference. Cambridge Univ Press (2000).

41

R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical
Computing, Vienna, Austria (2015).
Rosenbaum, P. “Interference between units in randomized experiments.” Journal of the American
Statistical Association, 102(477):191–200 (2007).
Rosenquist, J. N., Murabito, J., Fowler, J. H., and Christakis, N. A. “The spread of alcohol consumption
behavior in a large social network.” Annals of Internal Medicine, 152(7):426–433 (2010).
Ross, N. F. “Fundamentals of Stein’s method.” Probability Surveys, 8:210–293 (2011).
Rubin, D. “Comment: Neyman (1923) and causal inference in experiments and observational studies.”
Statistical Science, 5(4):472–480 (1990).
Shalizi, C. and Thomas, A. “Homophily and contagion are generically confounded in observational
social network studies.” Sociological Methods & Research, 40(2):211–239 (2011).
Shalizi, C. R. “Comment on “Why and When ’Flawed’ Social Network Analyses Still Yield Valid Tests
of no Contagion”.” Statistics, Politics, and Policy, 5(1) (2012).
Shalizi, C. R. and Rinaldo, A. “Consistency under Sampling of Exponential Random Graph Models.”
Annals of Statistics, 41(2):508–535 (2013).
Sobel, M. “What Do Randomized Studies of Housing Mobility Demonstrate?” Journal of the American
Statistical Association, 101(476):1398–1407 (2006).
Sofrygin, O. and van der Laan, M. J. “Semi-Parametric Estimation and Inference for the Mean Outcome
of the Single Time-Point Intervention in a Causally Connected Population.” U.C. Berkeley Division
of Biostatistics Working Paper Series, (Working Paper 344) (2015).
Sofrygin, O. and van der Laan, M. J. tmlenet: Targeted Maximum Likelihood Estimation for Network
Data (2015). R package version 0.1.0.
Sofrygin, O., van der Laan, M. J., and Neugebauer, R. simcausal: Simulating Longitudinal Data with
Causal Inference Applications (2015). R package version 0.5.0.
URL http://CRAN.R-project.org/package=simcausal
Stein, C. “A bound for the error in the normal approximation to the distribution of a sum of dependent
random variables.” In Proc. Sixth Berkeley Symp. Math. Stat. Prob., 583–602 (1972).
Tchetgen Tchetgen, E. J. and VanderWeele, T. “On causal inference in the presence of interference.”
Statistical Methods in Medical Research, 21(1):55–75 (2012).
Thomas, A. C. “The social contagion hypothesis: Comment on ‘Social contagion theory: Examining
dynamic social networks and human behavior’.” Statistical in Medicine, 32(4):581–590 (2013).
Tsiatis, A. Semiparametric theory and missing data. Springer Science & Business Media (2007).
Vallender, S. “Calculation of the Wasserstein distance between probability distributions on the line.”
Theory of Probability &amp; Its Applications, 18(4):784–786 (1974).
42

van der Laan, M. J. “Causal Inference for a Population of Causally Connected Units.” Journal of
Causal Inference, 0(0):2193–3677 (2014).
Van der Laan, M. J. and Robins, J. M. Unified methods for censored longitudinal data and causality.
Springer Science &amp; Business Media (2003).
Van der Laan, M. J. and Rose, S. Targeted learning: causal inference for observational and experimental
data. Springer Science & Business Media (2011).
van der Laan Mark, J., Polley Eric, C., et al. “Super learner.” Statistical Applications in Genetics and
Molecular Biology, 6(1):1–23 (2007).
Van Der Vaart, A. W. and Wellner, J. A. “Weak Convergence.” In Weak Convergence and Empirical
Processes, 16–28. Springer (1996).
VanderWeele, T. “Direct and indirect effects for neighborhood-based clustered and longitudinal data.”
Sociological Methods & Research, 38(4):515–544 (2010).
VanderWeele, T. and Tchetgen Tchetgen, E. “Bounding the Infectiousness Effect in Vaccine Trials.”
Epidemiology, 22(5):686 (2011a).
—. “Effect partitioning under interference in two-stage randomized vaccine trials.” Statistics & probability letters, 81(7):861–869 (2011b).
VanderWeele, T. J., Ogburn, E. L., and Tchetgen, E. J. T. “Why and when ’flawed’ social network
analyses still yield valid tests of no contagion.” Statistics, Politics, and Policy, 3(1) (2012).
Wasserman, S. “Comment on “Social contagion theory: Examining dynamic social networks and human
behavior” by Nicholas Christakis and James Fowler.” Statistics in Medicine, 32(4):578–580 (2013).
Watts, D. J. and Strogatz, S. H.

“Collective dynamics of ’small-world’ networks.”

Nature,

393(6684):440–442 (1998).
Young, J. G., Hernán, M. A., and Robins, J. M. “Identification, Estimation and Approximation of
Risk under Interventions that Depend on the Natural Value of Treatment Using Observational Data.”
Epidemiologic Methods, 3(1):1–19 (2014).

43

