Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
140,877 characters · 52 sections · 129 citation commands
Model-Based Inference and Experimental Design for Interference Using Partial Network Data
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 \fi
\if01 {
} \fi
\spacingset{1.9}
Interference occurs when one individual's treatment status impacts others' outcomes. Interference, also known as “spillover effects,” appears in multiple scientific domains, including the study of infectious diseases hudgens2008toward, tchetgen2012causal, studying peer influence manski1993identification, bramoulle2009identification, de2010identification, epple2011peer, goldsmith-pinkhami2013, public policy malani2021effect, imai2021causal, information diffusion banerjee2013diffusion, banerjee2019gossip, technology adoption beaman2021can, online platforms saveski2017detecting, pouget2018optimizing, pouget2019testing and online marketplaces ha2020counterfactual, johari2022experimental, among other domains.
Interference violates the stable unit treatment value assumption (SUTVA), which states that an individual's outcome is not impacted by the treatment status of their peers. When SUTVA is violated, each potential outcome, the counterfactual outcome under a given treatment assignment, could depend on all treatment assignments within the population. Valid inference for treatment effects under SUTVA violations is an active area of research, with solutions typically depending on a combination of exposure maps and structural causal models. Exposure maps categorize respondents according to their network characteristics and the vector of treatment statuses Aronow2017EstimatingExperiment, chandrasekhar2023general, while structural causal models identify specific pathways for influence between individuals van2012causal, ogburn2022causal.
Estimating causal effects under interference typically requires complete network data, which is expensive and onerous to collect or may not be available due to privacy constraints. Partially observed network data takes many forms: subgraph samples where a researcher observes the presence/absence for only a subset of possible connections, egocentric sampling using either specific links or aggregates, or network-based sampling methods such as snowball sampling or respondent-driven sampling. In each case, incomplete network information introduces miss-measurement in the exposure map. A person may have treated peers, for example, but if links to those peers are not observed, the researcher will think their outcome is totally orthogonal to the treatment.
This paper introduces a framework for estimation and inference of causal effects under partial network data arising from a single graph. Partial here means that we may observe some or no links or aggregate summaries of links, which we will formalize later. With such data, we recover multiple estimands including various conditional or average treatment effects. To do this, we define a broad class of structural causal models that are amenable to estimation using partial data. This class covers many empirically relevant schemas for interference, such as diffusion and its generalizations. Estimation leverages a dual approach: first, by using an iterated expectation method for de-biased estimation of model parameters with partial network data, and second, by managing the dependence of exogenous noise in the outcomes. chandrasekhar2011econometrics also introduced an iterated expectations strategy for cases where multiple networks are available and data are independent across networks. We tackle the more challenging inference task of single network asymptotics. Our method applies when the underlying graph has features captured by the class of node-exchangeable formation models, which we commonly see in practice and connect this to the problem of estimating effects of experiments. Previous methods breza2017using,breza2023consistently developed a related method to estimate network model features using a specific type of aggregated network data. Along with expanding to a wide range of partial data types, we extend this existing methodology to relax the requirement in previously published studies that traits be mutually distinct, a challenge to its usability in practice until now.
We also tackle the problem of experimental design associated with network exposure in scenarios where obtaining pristine network data ahead of randomized controlled trials (RCTs) is challenging or impossible. By collecting partial network data and employing a Bayesian optimization algorithm, we propose experimental designs that efficiently maximize treatment saturation tailored to specific estimands of interest. Our results demonstrate that this methodology not only surpasses traditional methods like inverse probability weighted (IPW) estimators in estimating global average treatment effects but also facilitates innovative seeding strategies that leverage the unique characteristics of partial network data. We demonstrate that these techniques can be used to assign treatment in such a way as to minimize estimator variance or to optimally seed for diffusion.
The remainder of the paper is structured as follows. We begin with a review of related work (section (ref). section (ref) defines the necessary background, then section (ref) describes the procedure for estimation and inference. section (ref) describes experimental design using partial network data and section (ref) provides empirical examples. We conclude in section (ref). Code to replicate the results in the paper is available at \url{https://github.com/SteveJWR/ardexp}, and an R package is available from \url{https://github.com/SteveJWR/SBMNetReg}.
We first provide a brief review of literature related to inference with partial network data, then move to an overview of causal inference under interference. Complete network data collection can be prohibitively expensive and restricted by privacy concerns breza2017using. Researchers typically work with partial network data derived from various sources such as survey samples, coarse geographic data, kinship information from censuses, or aggregated financial transactions. Comprehensive reviews of methods for handling network data can be found in de2017econometrics or graham2020network, and discussions on identification in network and related models are provided in manski2009identification. A direct approach is node subsampling, selecting a portion of nodes from the population and mapping the entire graph among them. If random sampling of nodes is infeasible, or if populations are sensitive or stigmatized, techniques like snowball or respondent-driven sampling offer a limited but focused view of the graph heckathorn1997respondent, goel2009respondent, goel2010assessing, baraff2016estimating, green2020consistency
When complete edge enumeration among node subsets is impractical, researchers adopt standard survey methods such as Aggregated Relational Data (ARD) collection. The main intuition is that each of the partial network designs mentioned above can be used to estimate a breakdown of each respondent's network in terms of observable characteristics. In ARD surveys, respondents are asked, “How many people do you know with trait X?" for various traits. Additional conditions may be added in addition to collect the type of connection that is relevant feehan2016quantity.Originally designed to estimate hard-to-reach populations like HIV-positive men in the US killworth1998estimation, scutelniciuc2012network, jing2014estimating, has been extended to a variety of other settings such as financial contagion models acemoglu2015systemic as well as more general network scale up methods utilized (NSUM) killworth1998social,kadushin2006scale,feehan2016generalizing,mccormick2020network and is notably 70 to 80% less costly than full network data collection breza2017using. Another standard survey method, egocentric sampling, asks respondents to consider specific individuals in their networks and provide detailed information about them, unlike the aggregate focus of ARD and is commonly used in applications such as contact tracing potter2011estimating, violence perpetration bond2017contagious or adolescent substance measurement huang2014peer.
The first task in causal inference problems, particularly in the presence of interference, is defining the target estimand. The global average treatment effect (GATE), for example, assesses the impact of treating everyone versus treating no one, considering peer effects Ugander2013GraphUniverses. Other interests might include the effect of specific treatment allocations, like identifying influential individuals kempe2003maximizing, banerjee2019gossip, often limited by policy constraints (e.g., subsidies for the ultra-poor as in anderson2007agricultural) or due to non-monotone peer effects dynamics banerjee2018less: treating everyone may change interaction dynamics in equilibrium. More generally Aronow2017EstimatingExperiment compare average treatment effects between two exposure configurations. A distinct but related line of work seeks to detect whether interference is present at all athey2018exact.
Models for peer influence like contagion jackson2008social, banerjee2013diffusion, beaman2021can, xiaoqi2023measuring or hearing models banerjee2019gossip structure interference analysis by identifying specific mechanisms that describe how connections between peers impact outcomes. auerbach2021local explore these effects through structural causal models focusing on nonparametric estimation, while our work emphasizes estimation, inference, and design using partial network data. Much of this literature assumes a fully observed graph, though a recent line of literature address imperfect or incompletely sampled graphs under certain conditions and for specific average causal effects Hardy2019EstimatingNetworks, yu2022estimating, cortez2022staggered. A related line of work examines sensitivity analysis for standard causal estimators under hidden treatment diffusion tortu2021causal.
Let $i \in \{1,2,\dots, n\} = \mathscr{V}$ denote a population of interacting individuals and let $\mathscr{G} = \mathscr{V} \times \mathscr{E}$ be the network by which interference is propagated; where $\mathscr{V}$ is the set of node vertices and $\mathscr{E} \subset \mathscr{V} \times \mathscr{V}$ is a set of edges (either directed or undirected). We can also extend this to weighted graphs, however binary networks are presented for simplicity. We can represent this graph by the adjacency matrix $G \in \{0,1\}^{n \times n}$. We consider binary treatments denoted by a treatment vector $\mathbf{a} \in \{0,1\}^n$ and let denote the potential outcome $Y_i(\mathbf{a}) \in \mathbb{R}$, under a treatment assignment $\mathbf{a}$, and $Y_i$ denote the actual observed outcome. Lastly, we assume that we have access to pre-treatment node-level covariates $X_i \in \mathbb{R} ^{m}$. In the remainder of the paper let $O$ and $o$ denote the usual big and little oh notation and $O_P$ and $o_P$ denote the stochastically bounded and convergence to $0$ in probability for sequences of random variables. We use $\tilde{O}$ if we are suppressing logarithmic factors in the rate. Let $||{\cdot}||_{p}$ denote and $p$-norm, and let $||{\cdot}||_F$ denote the Frobenius norm.
We use the framework of structural causal models, a nonparametric extension of structural equation models pearl2009causality. Similar approaches have been studied by ogburn2022causal and auerbach2021local in the case of fully observed networks. We derive a model that is amenable to estimation with partial data.
Let $Y_i(\mathbf{a})$ denote the potential outcome of $Y_i$ under a treatment allocation $\mathbf{a}$. The exposure mapping $V_i$ is represented as a function $f_V$ such that $V_i = f_{V}(\mathbf{a}, \varphi_i(G)) \in \mathbb{R} ^{p_V}$ where $\varphi_i$ is the relevant graph information for individual $i$ relative to their position with respect to treated individuals. We also allow for the potential outcome to be modulated by some additional confounder $S_i = f_{S}(\mathbf{X}, \vartheta_i(G)) \in \mathbb{R} ^{p_S}$. We model the potential outcomes $Y_i$ as a function of the exposure, type-value $S_i$ and some additional noise $\mathbf{\varepsilon}_{Y}$
The benefits of structural causal models are that they allow for the characterization of all causal effects in a system, as well as the distributions of counterfactuals. However, they require correct specification of the causal process, i.e. correct specification of the exposure map and the relevant confounders. Even if one can propose a model for interference, estimation is not straightforward due to the fact that we only observe partial graph information in $G^*$. Many common models of interference can be expressed as structural causal models, and can be thought of as parameterizations of $f_Y(S_i, V_i, \mathbf{\varepsilon}_{Y}) = f_Y(S_i, V_i, \mathbf{\varepsilon}_{Y}; \beta_0)$. This then reduces the challenge to estimating $\beta_0$ using partially observed data. The exogenous noise, $\epsilon_Y$, within our model is likely influenced by the graph structure, as interactions and peer effects can induce correlations in outcomes that extend beyond individual exposures. This complexity suggests that the noise, even if initially considered as external to the model, is intertwined with the network dynamics, reflecting the propagation and interference effects inherent in our structural causal framework.
We distinguish two types of target parameters. The first are the outcome model parameters, which parameterize the distribution of the outcome, exposure, and confounder $(Y,S,V)$. Specifically, $f_Y(S_i, V_i, \epsilon_Y) = f_Y(S_i, V_i, \epsilon_Y;\beta_0)$ under parameterization $\beta \in \mathbb{R}^p$. The true model parameters are $\beta_0 \in \mathbb{R} ^p$, identifiable through a moment equation $m$, $\mathbb{E}[m(Y_i,S_i,V_i, \beta_0)] = 0$, or through regression. In a simple diffusion model, this is the probability of infecting a neighboring node $q \in [0,1]$.
The second set of parameters we consider are the causal parameters, those involving the distributions of the counterfactuals. The main causal parameter we will consider is the expected average potential outcome on the complete network $G$, $\Psi(\mathbf{a}|G) = \frac{1}{n}\sum_{i = 1}^n\mathbb{E}[Y_i(\mathbf{a})]$, though these can also be made conditional on a covariate $x$: $\Psi(\mathbf{a}|x, G)$. Leveraging the structural causal model, we can define the these causal effects in terms of the structural causal model. We illustrate conditions for identification of these causal effects in section (ref). While we focus on defining causal quantities through conditional means, the nonparametric identification can also apply to other functionals like quantiles.
Inference for learning the causal parameters under the above assumptions now amounts to learning the distributional relationship between $Y_i$ and $S_i,V_i$. We consider settings where the assignment of treatments can be manipulated by an experimenter, which we discuss in section (ref). If one leverages this model, either through assumption or estimation, then we can use a structural causal model to generate expected potential outcomes under different treatment assignments $f_Y(S_i, V_i, \mathbf{\varepsilon}_{Y})$, which is precisely what is done in the case of seeding. A contrast of these frameworks is included in the Appendix in section (ref). The applicability of a model to a new population parallels challenges in distribution shift, as explored in shimodaira2000improving or wilkins2024multiply.
Adding structure to the potential outcomes model is standard in fields like economics, where researchers often propose models to explain how information or behaviors spread across networks. Many of these models include a temporal element. In our setting, we consider outcomes at a fixed time $T$, i.e., $Y_i(\mathbf{a}) = Y_{i,T}(\mathbf{a})$. For instance, banerjee2013diffusion explore a latent diffusion process in micro-lending, banerjee2019gossip study a hearing model for information diffusion, and beaman2021can analyze behavior adoption in agriculture through complex contagion. Additionally, centola2007complex differentiate the spread of information, often through single links, from behaviors that require multiple neighbors for network propagation.
A foundational model of information diffusion is based on simple contagion, and generalizations of SIR (Susceptible-Infected-Recovered) models on networks kermack1927contribution, giles1977mathematical. These models have been further studied and extended in various settings jackson2006diffusion, aral2009distinguishing, romero2011differences, chierichetti2011reconstructing, banerjee2013diffusion, banerjee2019gossip. Here we illustrate how the base model, under which many extensions are built, can be interpreted as a structural causal model. This interpretation can also be applied to complex contagion settings centola2007complex, beaman2021can, cencetti2023distinguishing.
Consider a scenario where initially infected (treated) seeds $\mathbf{a}$ transmit the infection to each neighbor with probability $q$ at each time-step $t \in {1,2,\dots, T}$, after which they are no longer infectious. An infection status at time $t$ is denoted as $Y_{it} = 1$. The overall outcome $Y_i= 1$ indicates whether a node was infected at any time up to $T$. For a simple case with $T = 2$, we model the transmission using Bernoulli random variables $\epsilon_{ij} \sim \text{Bernoulli}(q)$, representing potential infection from node $i$ to node $j$. Let $\mathbf{E}_{ij} = \epsilon_{ij}$ and $\mathbf{D} = \mathbf{E} \odot G$. This setup is depicted in Figure (ref). Given a random sample of the directed graph $\mathbf{D}$, one can characterize what would have happened if a node were treated, which is precisely the counterfactual. For instance in Figure (ref) we seed the left most which proceeds to propagate in steps $1$ and $2$. Additionally, one can construct the relevant exposure map for any fixed number of time steps $T$.
We consider several examples of exposure maps, though this list is not exhaustive.
Now that we've established our framework for interference, we next return to the data used for estimation. In our setting, we do not have access to the full graph $G$, but rather, have access to some summarizing function the graph $G^* = \zeta(G)$. tsiatis2006semiparametric uses the term coarsened data to refer to such partial measurements of missing data in general, not necessarily in the network setting. “Coarsened” is apt in our setting because our method works on partial graph structure that give (either directly or by construction) an estimate of linking rates across population members with different combinations of observable traits. A non-exhaustive set of examples of partially measured network data include induced subgraphs or egocentric sampling freeman1982centered, almquist2012random, respondent driven sampling heckathorn1997respondent, aggregated relational data killworth1998estimation, respondent driven sampling heckathorn1997respondent, goel2009respondent, goel2010assessing, green2020consistency and more.
In order to infer about the distribution of the missing part of the graph, we propose that $G \sim \theta_0$ where we assume that $\theta_0 \in \Theta$ denotes the parameters of a random graph model. In this case, for each $i$, there is an a latent $\xi_i$ parameter such that $P(G_{ij} = 1 |\xi_i, \xi_j) = \tilde g(\xi_i, \xi_j)$ for some function symmetric, measurable, $\tilde g$, known as a graphon lovasz2006limits, orbanz2015bayesian. Many common graph models, such as latent space models Hoff2002LatentAnalysis, handcock2007model, lubold2023identifying, wilkins2022asymptotically, are included in this category. Graphons are appealing in this context because, following airoldi2013stochastic and gao2015rate, they can be approximated arbitrarily well using latent types assigned to each node. Said another way, graphons introduce complex dependence in the network-generating mechanism through clustering induced by latent types associated with each node. In our inferential procedures in section (ref), the general procedures involve estimation from a missing data perspective. This will involve estimating the graph model $\hat \theta := \hat \theta(G^*)$ then inferring about the distribution $G|G^*, \hat \theta$. Further details for estimating the graph model are included in section (ref).
We first illustrate an identification procedure for the causal effect without a-priori imposing any model structure. These are analogous to standard causal identification assumptions, adapted to our framework.
Conditioning the graph confounder $S_i$ captures the heterogeneity of the outcomes when observed with a given exposure. In simple contagion models, nodes are equivalent, and this independence occurs naturally without conditioning. section (ref) discusses an example from ugander2023randomized where conditioning on node degree suffices for any randomization.
This assumption can be simply understood as the exposure is correctly specified.
This assumption states that once we have adjusted for $V_i$ and $S_i$, then the potential outcomes are independent of the network $G$. These assumptions allow us to express the causal estimand through observational data.
For brevity, we denote the true conditional mean $\mathbb{E}[Y_i|V_i = v, S_i = s]$ as $h_0(s,v)$ and denote $h(s,v;\beta)$ a model to estimate $h_0(s,v)$. Given a network model $\theta$, observed graph data $G^*$, and a conditional model $h(s,v;\beta)$ we can also define the expected average treatment effect
where under the correct model conditional model and graph model $\mathbb{E}[\Psi(\mathbf{a}| G)|\mathbf{a}, \mathbf{X}, \theta_0] = \Psi(\mathbf{a}| \beta_0, G^*, \theta_0)$. In Appendix (ref), we illustrate when this population average effect under any draw of the network $\Psi(\mathbf{a}| G)$ will be close to the average over the model class $\Psi(\mathbf{a}| \beta_0, G^*, \theta_0)$; and study plug-in estimators $\hat \Psi(\mathbf{a}| G) = \Psi(\mathbf{a}| \hat \beta, G^*, \hat \theta)$.
We outline our method for estimating parameters with partial network data. Developing these results requires two theoretical tools: a fast estimation rate for network model parameters $\theta_0$, and a suitable central limit theorem for scenarios with correlated outcomes. Outcome Model Parameters and Estimators Next we consider estimating the outcome model parameters $\beta_0$. We present two methods for estimating such parameters, instrumentation in a linear model, and $Z$ estimators. The iterated expectation procedure for estimating such parameters was introduced in chandrasekhar2011econometrics, however, we extend inference to the single network setting. Similar approaches exist for peer effects models boucher2020estimating.
We first illustrate identification of the conditional model under a linear model assumption.
where $\mathbb{E}[\varepsilon_i] = 0$ and there can be general correlation $\mathrm{Var}[\mathbf{\varepsilon}] = \Sigma$. Without access to the network data, one can recover the model parameters through conditional expectation
where we create a new set of features $\tilde H_i= \mathbb{E}[\tilde h(S_i(G), V_i(\mathbf{a}, G))| \mathbf{a}, \mathbf{X}, G^*, \theta_0]$ by averaging over the network model. Here identification comes from the variation of these averaged features $\tilde H_i$ over the population. More clearly, letting $\mathbf{\tilde H} \in \mathbb{R} ^{n \times p}$ denote the design matrix of this model, identification comes from the linear independence of the columns of $\mathbf{\tilde H}$.
In other cases, parameters may be defined through a moment equation, and can be used to construct a $Z$-estimator, for example, generalized linear models (GLMs). These parameters can be identified using an estimating equation approach where given a moment function $\tilde m(Y_i,S_i,V_i; \beta)$ such that $\mathbb{E}[\tilde m(Y_i,S_i,V_i; \beta)|\mathbf{a}, \mathbf{X}, G] = 0$ if and only if $\beta = \beta_0$. Through the use of iterated expectations, we can define a new estimating equation, by marginalizing over the draws of the graph model then applying iterated expectations
Identification arises from the variation of exposure and confounders, such that $\beta = \beta_0$ if and only if $\mathbb{E}\left[\frac{1}{n}\sum_{i = 1}^n m_i(Y_i, S_i, V_i; \beta, \theta_0)\bigg| G^*, \mathbf{a}, \mathbf{X}, \theta_0\right] = 0$. Exact conditions vary by parameter, but GLMs can use a similar strategy as linear models.
We introduce a general procedure for estimating the outcome model parameters. We also illustrate inference for estimation of a causal target parameter on a particular graph $G$. We present a pseudo-code approach to the procedure in Algorithm (ref). Let $\tilde Z_i= (Y_i, S_i, V_i)$ denote the full (including unobserved) data, and let $\mathbf{Z} = (\mathbf{Y}, \mathbf{a}, \mathbf{X}, G^*)$ denote the observed data.
Step (ref) asks the practitioner to propose a response model given the treatment, i.e. the causal model in section (ref). Step (ref) estimates the generative model given the partial network data and the node covariates observed. We give theoretical results where the formation model is a stochastic blockmodel, then give rate estimation relative to the more general graphon approach. Step (ref) estimates the parameter by marginalizing the estimating function over the graph model. Lastly, Step (ref) is optional if the target parameter is a plug-in estimator of the causal parameter using the regression model. We discuss inference for the plug-in estimate of causal parameters using a delta method argument in the Appendix (ref). We next give our asymptotic results, then provide an example of this algorithm in section (ref).
The asymptotic results for both the Z-estimator and the linear model will depend on being able to establish a central limit theorem based on the exogenous noise. To establish asymptotic properties for our outcomes on a network, we extend the application of the central limit theorem (CLT) to structures not commonly associated with traditional time series or spatial dependencies. Nonetheless when the exogenous noise is correlated, we will need a method of handling the central limit theorem. Specifically, we utilize a general version of the CLT for dependent data from chandrasekhar2023general. For brevity in presentation, we leave the full detail of this central limit theorem to the appendix.
We denote $g_i(\mathbf{Z}; \beta) = m_i(\mathbf{Y}; \mathbf{a}, \mathbf{X}, \beta, G^*,\theta_0)$ to be the moment function evaluated using the true generative model and correspondingly $g_n(\mathbf{Z}; \beta) = \frac{1}{n}\sum_{i = 1}^n g_i(\mathbf{Z}; \beta)$. Further, define the (normalized) random vector of the estimating function evaluated at the correct model parameters $\mathcal{E}_i = \frac{1}{n}g_i(\mathbf{Z}; \beta_0)$. And lastly let $D_n(\mathbf{Z}; \beta_0) = \nabla_{\beta} g_n(\mathbf{Z}; \beta_0) \in \mathbb{R} ^{p \times p}$ denote the gradient of the estimating equation $g_n(\mathbf{Z}; \beta)$. To develop valid inference, we must estimate the graph model quickly enough to disregard the graph estimation component during inference. We will next present the theorem and discuss the assumptions further.
The first set of assumptions ensures the consistency of $Z$-estimators, typically derived from uniform laws of large numbers as discussed in andrews1987consistency,newey1994large. The second set involves conditions that make the graph model's estimation negligible, requiring the estimating functions to be smooth with respect to the graph parameters.
The final set of assumptions, stated in (ref), are utilized so that $\mathcal{E}_i$ satisfy a central limit theorem chandrasekhar2023general. This assumption is required if the data exhibit further dependence after controlling for graph parameters (if, for example, there are latent factors that impact both outcomes and the propensity to form ties). The main idea of chandrasekhar2023general is to represent dependence in terms of “affinity sets” where the majority of dependence structure is captured within sets, leaving little between sets. In the modelling of social behaviours beyond just considering outcomes as a function of the exposure observed, outcomes may be further correlated, beyond examples of spatial dependence or heteroskedasticity. In practice we can include these dependencies through correlation terms matching the generative graph model, such as between blocks of a stochastic blockmodel or via latent positions in a latent space model.
Here, $r(n)$ describes the effective rate at which the variance converges. For the estimation of the graph model $\theta_0$ to be considered negligible, it must occur more rapidly than $r(n)$. In cases of independent or minimally dependent noise, it is typical for $r(n) \approx n^{-1/2}$. Alternatively, in different scenarios, $\mathcal{E}_i$ might exhibit correlation within densely connected blocks of the network, such as during a diffusion process in a stochastic blockmodel with $k_n$ densely linked blocks (refer to chandrasekhar2023general, section 4.4). In such cases, $r(n)$ is generally on the order of $k_n^{-1/2}$. If both $r(n)$ and $s(n)$ approach zero, but $\frac{s(n)}{r(n)}$ diverges or stabilizes at a nonzero constant, a consistent estimator for the outcome model parameters can still be obtained. However, its asymptotic distribution may be influenced by the graph model estimation, necessitating a tailored inference approach.
An analogous argument follows when conducting inference using a linear model. For the sake of brevity and avoiding repetition, we include it in the Appendix in section (ref). In Theorem (ref) we present a summary.
We next discuss the estimation of the generative model for the network using a variety of data types. We first demonstrate results for estimating parameters in a stochastic blockmodel. We then extend the result to view the blockmodel as an approximation of a graphon. breza2017using and breza2023consistently consistently estimate a generative model for ARD with mutually exclusive traits. We extend this work by introducing a novel method for estimating the stochastic blockmodel with non-mutually exclusive traits using constrained least squares approach. This innovation has major implications for practice since it can dramatically reduce survey length (since asking about multiple categories requires constructing separate questions for each intersection to make the traits mututally exclusive). Our approach applies to a range of partial network data, not just ARD, but we summarize the resulting rates using ARD for a variety of model classes in the appendix in Table (ref) for comparison with existing literature.
In the main text, we concentrate on estimating the stochastic blockmodel using ARD. In the Appendix in section (ref) we estimate generative models using partial network data such as subgraph sampling and develop similar rates for the stochastic blockmodel for subgraph sampling and reference a similar result for respondent driven sampling.
Recall that $X^*_{it}$ represents a set of ARD response vectors. breza2023consistently show that we can consistently estimate the connection probabilities between latent types, however, we present an improved version of the SBM estimator which allows for an non-mutually exclusive traits. Let $n_t$ denote the total number of individuals of trait type $t$. Let $N'_k$ denote the nodes in our sample in group $k$, and let $n_k$ denote the number of nodes in the graph in group $k$. We cluster the node memberships according to Algorithm (ref).
After we obtain a clustering, we can estimate the stochastic blockmodel. Let $\hat \Omega_{kt} = \hat N_{kt}/N_t$ where $N_{kt}$ are the number of traits in the estimated group $k$ and with trait $t$, and $N_t$ are the number of individuals with trait $t$, and $\Omega_{kt} = N_{kt}/N_t$, the analogous population quantity. We next define the probability matrix of observing a connection of group $k$ with a trait $t$. $\mathbf{\tilde P}_{kt} = \sum_{k'} \mathbf{P}_{kk'}\omega_{k't}$, where $\mathbf{\tilde P}_{kt} = P(G_{ij} = 1|k_j = k, t_i = t)$. This relationship can be expressed in a linear system $\mathbf{\tilde P} = \Omega \mathbf{P}$ where $\Omega \in \mathbb{R}^{T \times K}$ and $\Omega_{kt} = \omega_{kt}$. If $\Omega$ is of full column rank, then a unique solution will exist as:
In general, one can symmetrize $\mathbf{\hat P}_{kk'}$ after the estimate to ensure the constraints of an undirected stochastic blockmodel are satisfied. Alternatively, once can also minimize the constrained least squares objective which can be implemented using standard convex solvers such as CVX Fu2020CVXR:Optimization
breza2023consistently develop a method for consistently estimating the stochastic blockmodel. We extend their result by obtaining a rate for estimating model parameters (Lemma (ref)) and relax the assumption that of mutually exclusive traits. We differentiate between the estimated cross-group probabilities $\mathbf{P}^{(\mathbf{\hat k})}$ and those under known membership $\mathbf{P}^{(\mathbf{k})}$.
We contrast our results to the optimal estimation rate for a stochastic blockmodel from gao2015rate, $\tilde O_P(n^{-1/2})$. Our rate appears faster due to the complexity difference in clustering problems. Our clustering benefits from node-level traits, which provide extra information. As the network grows, the normalized ARD vector converges to its mean, simplifying clustering and resulting in a faster rate.
We use a stochastic blockmodel as it effectively approximates a general graphon class. Even if $\theta_0$ belongs to a smooth graphon class rather than a stochastic blockmodel, we can still bound the bias in estimating the relevant model parameters. Consider a scenario where edges are generated under a true graphon model $\tilde g$ where $\eta_{ij} = \tilde g(\xi_i, \xi_j) = P(G_{ij} = 1|\mathbf{\xi}) \text{ where } \mathbf{\xi} \sim_{iid} P_{\mathbf{\xi}} \in [0,1]. $ Let $\mathcal{H}_{\alpha}(M)$ denote a smooth graphon class defined via the $\alpha$-$M$-H\"{o}lder class as follows. Let $\mathcal{D} = [0,1]^2 \cap x \leq y$ denote the domain of $(x,y)$. We define the norm $||{\tilde g}||_{\mathcal{H}_\alpha}$ as: $$ ||{\tilde g}||_{\mathcal{H}_\alpha} = \max_{j + k \leq \lfloor \alpha \rfloor} \sup_{x,y \in \mathcal{D}} |\nabla_{jk} \tilde g(x,y)| + \max_{j + k = \lfloor \alpha \rfloor}\sup_{(x,y) \not = (x'y') \in \mathcal{D}}\frac{\nabla_{jk} \tilde g(x,y) - \nabla_{jk} \tilde g(x',y')}{(|x- x'| + |y - y'|)^{\alpha - \lfloor \alpha \rfloor}}$$ and the H\"{o}lder class corresponding to this norm as $$ \mathcal{H}_\alpha(M) = \{||{\tilde g}||_{\mathcal{H}_\alpha} \leq M: \tilde g(x,y) = \tilde g(y,x); 0 \leq \tilde g(x,y) \leq 1\}. $$ Prior work has focused on the approximability of a stochastic blockmodel to any element of a smooth graphon class. In particular there will always be some assignment of block memberships such that we can bound the $2$-norm probability deviation from the true model.
In practice, we don't directly select clusters; misspecified clusters relate to observed traits, thus this bound holds only under good alignment of clusters. This bound is a worst-case scenario and may be overly conservative regarding observed bias. Future work could involve sensitivity analysis of the response function and the latent graph model.
So far, our focus has been on estimating model parameters given a treatment assignment $\mathbf{a}$. We now explore experimental design methods that leverage partial network data to choose $\mathbf{a}$ to minimize the variance of our estimands. Leveraging partial network data for this purpose is particularly appealing in practice, since it requires substantially less investment than collecting full network data and could be collected as part of creating a sampling frame in settings where researchers collect data to construct the frame.
We consider saturation randomization experiments, which divide the dataset into $J$ clusters of size $n_j$. A proportion $\tau_j$ of each cluster is assigned the treatment, totaling $n_t = \sum_{j = 1}^{J} \tau_j n_j$, and generally will not be the same “blocks" as those in a graph model if the graph model uses discrete factors (e.g. stochastic blockmodel). Practically, due to budget constraints, the set of possible saturation levels $\mathbf{\tau}$ is limited to $\mathcal{T} \subset [0,1]^{J}$. For example, this could be due to limited resources like a finite vouchers in a vaccine trial.
Our goal is to optimize the asymptotic variance of a function of the model parameter $\hat \beta$ in section (ref). We highlight this by optimizing the variance of the estimates of linear contrasts of the parameters $\phi^T \beta$. When using the stochastic block model for the network model these treatment blocks could align with the model blocks, however this need not (and likely won't) be the case. They could, instead, be based on observed characteristics (e.g. geography, classrooms).
Denote the variance of the target contrast parameter conditional on the treatment assignment as $\mathbf{a}$: $\upsilon^{\phi}(\mathbf{a}; \theta) = \mathrm{Var}(\phi^T \hat \beta|\mathbf{a}, \theta)$. Ideally, the goal is to find a treatment assignment $a^*$ that minimizes the variance of the contrast: $a^* = \operatorname*{arg\,min}_{\mathbf{a} \in \{0,1\}^n} \upsilon^{\phi}(\mathbf{a}; \theta).$ Without added structure, optimizing treatment assignments is NP-hard, requiring a search over $2^n$ possible assignments. By changing the objective to one where we optimize over a set of saturation levels over a set of groups $\mathbf{\tau} \in [0, 1]^{J}$, we simplify the problem so that it is no longer NP-hard (i.e. since $J \ll n$ typically) and is therefore tractable. The distribution of treatment assignments, $\mathbf{a}$, under $\mathbf{\tau}$ is denoted by $P_{\mathbf{\tau}}$, and we aim to minimize: $$\mathcal{V}(\mathbf{\tau}; \theta) = \mathbb{E}_{\mathbf{a} \sim P_{\mathbf{\tau}}}[\upsilon^{\phi}(\mathbf{a}; \theta_0)].$$
In Algorithm (ref), we present a method for evaluating the variance of a linear model using a generic feature map $\tilde h$ for a given treatment assignment $\mathbf{a}$ and a graph model $\theta$. A general approach for Z-estimators is detailed in the Appendix. Algorithm (ref) operates under specific assumptions about the covariance matrix $\Sigma$, which may include correlations within densely connected network components. We will present our algorithm for minimizing this variance using Bayesian optimization, which accounts for the uncertainty in the outcome, given a graph model. In the appendix we give an extension which also which incorporates network model uncertainty $\hat \theta$ (section (ref)).
Bayesian Optimization. Calculating the average variance $\mathcal{V}(\mathbf{\tau}; \hat \theta)$ in Algorithm (ref) is computationally intensive to evaluate and often non-convex. Since the number of cluster saturation tends to be relatively small, this suggests that Bayesian optimization is an appropriate method for minimizing this saturation variance. Let $\mathcal{V}(\tau) := \mathcal{V}(\tau; \hat \theta)$ denote our objective function of the variance evaluated using an estimate of the network model $\hat \theta$. Given a set of pilot points $\mathbf{\tau}_1, \mathbf{\tau}_2, \dots, \mathbf{\tau}_{n_0}$ (i.e. uniformly sampled on $\mathcal{T}$) we propose a Gaussian process prior satisfying
where $\mathrm{Cov}[\mathcal{V}(\mathbf{\tau}_{i}), \mathcal{V}(\mathbf{\tau}_j)] = \Sigma_0(\mathbf{\tau}_i,\mathbf{\tau}_j)$ where $\Sigma_0$ is a positive semidefinite kernel function. As a default, we use the Gaussian kernel $\Sigma_0(x,x') = \alpha_0 \exp\left( -||{x - x'}||^2\right)$. We can then use this prior to define a posterior over remainder of the design space $\mathcal{T}$
From this posterior, we define an acquisition function $A(\mathbf{\tau})$. As a default, we choose the upper confidence bound (UCB) acquisition function $A(\mathbf{\tau}) = \mu_n(\mathbf{\tau}) - \kappa \sigma_n(\mathbf{\tau})$ for a chosen $\kappa$ (where we set $\kappa = 2$). This method is implemented in the R package rBayesianOptimization, which uses GPfit rBayesianOptimization_package, macdonald2015gpfit. For a detailed review of Bayesian optimization techniques, refer to frazier2018tutorial. We evaluate the complete Bayesian optimization procedure in Algorithm (ref), where we apply the procedure for $N_0$ iterations.
The quality of optimization over $N_0$ iterations depends on the smoothness of $\mathcal{V}(\mathbf{\tau})$. Since variance might diverge under some settings (e.g., as $\mathbf{\tau} \to 0$), a simple alternative is to maximize $\exp(-\mathcal{V}(\mathbf{\tau}))$ instead. The closeness of the maximizer after $N_0$ iterations hinges on the smoothness of $\exp(-\mathcal{V}(\mathbf{\tau}))$, which we assume belongs to a reproducing kernel Hilbert space, $\mathcal{H}$, with a bounded kernel $\Sigma_0(x,x') \leq B$. This function's smoothness affects the approximation rate, detailed in srinivas2009gaussian. For instance, with Gaussian kernel $\Sigma_0$, the approximation error is $\exp(-\mathcal{V}(\mathbf{\tau}^*)) \geq \frac{1}{N_0}\sum_{m = 1}^{N_0}\exp(-\mathcal{V}(\mathbf{\tau}_m)) + O_P(\frac{B\sqrt{\log(N_0)^{K + 1}} + \log(N_0)^{K + 1}}{\sqrt{N_0}})$. Similar findings apply to Matern and linear kernels per srinivas2009gaussian.
Given a model of the potential outcomes, we may also leverage this model for optimal seeding, a task that is NP-hard kempe2003maximizing in general. Many contagion models are exchangeable given an exposure, and with only block information available, then we can reduce our search space to that over block saturation. In our case, where exact network structures are unknown, we determine the optimal blocks for seeding. When $K \ll n$, this structure significantly reduces computational efforts, and we only need to decide how many seeds to allocate to each of the $K$ clusters.
The model leveraged for the outcome $f_Y(V_i, S_i, \mathbf{\varepsilon}_Y)$ could be a predefined model based on domain knowledge, such as complex contagion used by beaman2021can. In other scenarios, this might be estimated (e.g., simulation using $f_Y(V_i, S_i, \mathbf{\varepsilon}_Y; \widehat \beta)$ in place of $f_Y(V_i, S_i, \mathbf{\varepsilon}_Y)$). This is demonstrated in Algorithm (ref) (line 5).
When the total number of seeds (see Algorithm (ref), line 3) is small, it is computationally feasible to implement exactly. Alternatively, we could use Bayesian optimization to control treatment saturation levels.
In this section, we present three empirical examples to illustrate our framework's utility in estimating causal effects, designing experiments, and implementing seeding strategies. We adopt a semi-synthetic approach in our examples, where the outcomes are simulated based on processes derived from real networks. The networks analyzed pertain to observational and experimental studies focused on information diffusion in rural villages in India and Malawi, as discussed in banerjee2013diffusion, banerjee2019gossip, and beaman2021can. These networks consist of 30-400 households per village. To ensure continuity across the examples, we generate ARD as the partial data type and model the networks using stochastic blockmodels for each case, however the use of other network generative models and partial network datatype are applicable in these cases.
When covariates are available for all nodes, we use them to construct ARD. If covariates are missing, we apply the Leiden algorithm traag2019louvain in igraph csardi2006igraph to cluster the network and treat these clusters as traits. Table (ref) details which datasets used actual traits versus clustering to manage trait numbers in our simulations.
In Section (ref), we use networks from banerjee2013diffusion, which include various social relations from 70 villages, each with 80 to 350 individuals per household. In Section (ref), we use networks from banerjee2019gossip, consisting of 68 similar-sized villages, and repeat the simulations 500 times per village. In Section (ref), we include networks from beaman2021can, excluding those with insufficient connections for diffusion. This leaves 114 villages with 30 to 350 households per village, and we repeat the simulations 2000 times per village.
In our first example, the aim is to estimate the global average treatment. We consider the example from ugander2023randomized and generate a set of potential outcomes according to the following model $$ Y_i(\mathbf{0}) = \frac{d_i}{\bar d}\cdot \left(\alpha + bX_i + \sigma \epsilon_i \right), \quad Y_i(\mathbf{a}) = Y_i(\mathbf{0})\cdot\left(1 + \delta a_i + \gamma\frac{\sum_{j \in [n]} G_{ij}a_j}{d_i} \right) $$ where $\epsilon_i \sim_{iid} N(0,1)$ is some independent noise, and $X_i$ is a covariate that varies throughout the network, $d_i$ is the degree of individual $i$ and $\bar d$ is the average degree across the network. We set $\alpha = 1$, $b = 1$, $\delta = 1$, $\sigma = 0.5$ and $\gamma = -0.5$. The global average treatment effect in this model is $\frac{1}{n}\sum_{i = 1}^nY_i(\mathbf{0})(\delta + \gamma) = \Psi(\mathbf{a} = 1|G) - \Psi(\mathbf{a} = 0|G)$. The exposure is the individual treatment in conjunction with the average treatment of neighbors, and the graph confounder include the degree ratio and node level covariates $$f_V(\mathbf{a}; \varphi_i(G)) = \left(a_i, \frac{\sum_{j \in [n]} G_{ij}a_j}{\bar d}\right), \quad f_{S}(\mathbf{X};\vartheta_i(G)) = \left(\frac{d_i}{\bar d}, X_i\right). $$
We evaluate the effectiveness of graph cluster randomization by comparing a Horvitz-Thompson estimator Ugander2013GraphUniverses to a difference in means estimator under a cluster randomized design. In this design, half of the clusters receive no treatment (saturation of $0$) and the other half receive full treatment (saturation of $1$). We vary the number of clusters from 4 to 16 but display results only for 4 and 10 clusters in Figure (ref) for clarity.
Figure (ref) shows that the full data regression model performs the best, as it leverages more information than the ARD approaches. However, the ARD version still effectively minimizes bias (Figure (ref)) and RMSE (Figure (ref)). In our simulations of dense graphs with few clusters, the Horvitz-Thompson Estimator faces challenges as the network grows—almost all nodes have at least one neighbor with a treatment different than their own. The difference in means estimator shows consistent bias, due to not using heterogeneous covariate information. While regression with complete data is most effective, using partial network data still yields comparably good results.
We next highlight aspects of experimental design using an information diffusion example based on the hearing model referenced in section (ref). At each time step the previously infected nodes are susceptible again the nodes infected in the last round will infect their neighbors with probability $q_{t + 1}$. We repeat this for $T = 3$ rounds. Let $N_i$ denote the total number of infections after the process. We then sample some binary response $P(Y_i = 1|N_i) = \text{logit}(\alpha_0 + \alpha_1 N_i)$ where $\alpha_0$ and $\alpha_1$.
In this case, $V_i = \mathbb{E}[N_i|\mathbf{a}] = \sum_{t = 0}^3 \beta_t\mathbf{a}(G^{t})_{i}$ where $\beta_t = \prod_{j = 1}^t q_j$. We estimate the coefficients in each of these cases letting $V_i = \mathbb{E}[N_i|\mathbf{a}]$ be the exposure mapping. We then generate the outcomes according to the exposure received
where $\Lambda(\cdot)$ is the logistic function. For our experiments, we set $\beta = (0,0.5,0.05,0.005)$.
In the dataset, seeds are assigned uniformly with either 3 or 5 seeds per network. Following our procedure in section (ref), we compute the optimal seed allocations, ensuring no cluster receives more seeds than available in the actual experiment (either 3 or 5). In practice our Bayesian optimization procedure starts by randomly sampling the target space 20 times, followed by 20 iterations to refine saturation. We then compare the estimates for $\alpha_1$ and all model parameters as shown in Figure (ref). The results indicate that a more strategically designed experiment generally yields more significant gains than directly using the graph parameters. On average, using optimized designs rather than uniform random designs when collecting network data significantly reduced RMSE. Specifically, for estimating $\alpha_1$, the optimized design decreased RMSE by 38% ($\pm$12%) compared to 11% ($\pm$2%) with complete data (where the brackets refer to the 95% confidence interval of the mean estimate across simulations). For all parameters, the optimized design resulted in a 45% ($\pm$10%) reduction in RMSE, versus an 18% ($\pm$2%) reduction with complete data.
We apply our methodology to the seeding problem described in beaman2021can, where the diffusion of pit-planting technology among Malawian farmers follows a complex contagion process. The outcome model is defined as $Y_i = f_Y(S_i, V_i, \mathbf{\epsilon}_Y)$, with individuals having a threshold $\varsigma_i \sim N_{[0,\infty)}(\lambda, 0.1)$ for spreading infection based on neighbor infections from the previous time (where $N_{(a,b)}(\mu,\sigma)$ refers to the $\mu,\sigma$ normal distribution truncated on the interval $(a,b)$) . This process is simulated over three time periods to align with their experimental design, setting $\lambda = 2$ and repeating 2000 times for $K = 8$ clusters to determine optimal seeding groups.
We explore two seeding strategies: randomly assigning seeds to the top two members of optimal clusters, and seeding the nodes with the highest degrees within these clusters. We compare these strategies to common degree targeting, noting that our max degree method typically yields the highest adoption rates, especially in larger, sparser villages, as illustrated in Figure (ref). However, in very small or dense networks, the performance differences between strategies are negligible. Across all graphs we find the optimal seeding strategy to increase adoption by $1.50$ $(\pm 0.16)$ times relative to degree seeding, while the optimal blocks was $1.13$ $(\pm 0.12)$ times and optimal degree within blocks increased adoption by $1.28$ $\pm (0.13)$ times.
We introduce a framework that identifies causal effects under interference using a structural causal model, facilitating inference with partial network data. The framework is general and can be applied using broad class of outcome models and graph models. Our outcome modelling approach leveraging node-level heterogeneity and exposure mappings allow for the estimation of all causal effects, rather that other methods which tend to focus on a single causal effect like the GATE. Demonstrations through semi-synthetic problems highlight its effectiveness, matching or surpassing fully observed data methods in certain scenarios.
Our method highlights that directly modeling interference mechanisms offers several advantages, including leveraging transportability of outcome models for seeding and inference for experimental designs when estimating effects under interference.
Future studies might consider semiparametric approaches to estimation with partial data like those in auerbach2022identification. Additional structured assumptions on potential outcomes as suggested in belloni2022neighborhood could also be explored. Currently, our focus has been on analyzing problems at a single time point. However, future research could extend to designing experiments with panel data and staggered rollouts. It would also be worthwhile to develop classes of outcome models that more explicitly incorporate this temporal structure.
\numberwithin{equation}{section} \setcounter{section}{0}
We contrast the approaches of a fixed outcome approach as in Aronow2017EstimatingExperiment to a structural causal model approach. In the former approach, each individual has a distinct outcome under an exposure $v$, $Y_i(v)$. Though such an approach is robust for learning parameters such as average treatment effects $\frac{1}{n} \sum_{i = 1}^n Y_i(v)$, the information in an individual $i$'s potential outcome is completely distinct from individual $j$. This important details has important downstream implications.
Consider the simple contagion model from the example in section (ref) which takes place in a single time period ($T = 1$). Consider the nodes $i,j$ in Figure (ref) with seeded nodes in blue. Suppose that at time $T = 1$, that each neighbour of a treated node is infected with probability $q$. Since each one has only a single treated neighbor the distribution of the infection probability $P(Y_i = 1|\mathbf{a}, G)$ $i$ and $j$ are equivalent as their exposures are identical (i.e. they are each connected to a single seed node). However, in the finite sample framework the potential outcomes of any two nodes with a single treated neighbor can be arbitrarily different ($Y_i(v) \not = Y_j(v)$).
This nonparametric structure imposed on the potential outcomes later imposes restrictions on the degree of influence of others a node can have for estimation, thereby limiting this framework to examples with local dependencies (a phenomena also seen in ogburn2022causal).
In many nonparametric approaches to estimating causal quantities under interference, inverse probability weighted (IPW) estimates can be developed given a randomization scheme, i.e. distribution of the assignments $P(\mathbf{a})$ Aronow2017EstimatingExperiment. This is useful as it can be used to develop estimators for causal effects any exchangeability assumptions on the potential outcomes. However when $V_i$ is not observed directly, we must leverage additional structure in order to estimate any causal effects.
Our objective is to understand the model's structure and often apply it to tasks such as seeding. Thus, we rely on a correct model specification. The challenge with developing an IPW estimator arises when exposure is not observed. In such cases, it becomes impossible to determine which potential outcome was observed, violating the causal consistency assumption. Specifically, we don't know which potential outcome $Y_i$ represents (i.e., which exposure $v, Y_i = Y_i(v)$).
In this section we discuss extensions to several aspects of the paper with respect to the paper. Before proceeding we also introduce the full statement of the generalized central limit theorem result which we use to derive our asymptotic results.
Although network models do not neatly fit into conventional time series or spatial dependency categories, we provide a general framework by satisfying the necessary conditions through common dependence assumptions such as $M$-dependence. This includes scenarios characterized by $\alpha$-, $\phi$-, or $\rho$-mixing bradley2005basic. Our approach begins by defining affinity sets, (sets for which there is high correlation with an outcome) that form the foundational framework for applying the CLT, setting the stage for demonstrating its relevance and utility in analyzing network data.
The affinity sets can be used to construct a matrix which contains the bulk of the covariance across observations and dimensions. The regularity conditions can be understood as control of the covariance within affinity sets (ref), control of the covariance across affinity sets (ref) and control of the covariance outside of the affinity sets (ref). We collectively refer to these as the affinity set conditions. The affinity sets can be used to construct a covariance matrix $\Gamma_{n,dd'} = \sum_{i = 1}^{n}\sum_{(j,d') \in \mathcal{A}_{(i,d)}^{(n)}} \text{cov}(W^{(n)}_{i,d},W^{(n)}_{j,d'})$.
The authors illustrate several examples under which these conditions are sufficient for the this central limit theorem to hold. This theorem will be useful for proving our asymptotic results.
Here we first give the full theorem and regularity conditions with respect to the linear model.
We next discuss the estimation of generative models of network formation using several datatypes. We summarize the information for using ARD in Table (ref) as discuss similar rates for other datatypes.
We illustrate that it is possible to estimate the stochastic blockmodel using a diverse set of partial and sampled network data types. In each case, $\mathbf{P}_{kk'}$ refer to the cross-block probabilities, while $k_i \in \{1,2,\dots, K\}$ denote the node memberships. We consider partial network data to be any subset of the network data which can be used to generate an estimate of the generative model $\widehat \theta$.
Lastly, we discuss respondent driven sampling. In this setting, community membership can be defined based on a partition of the covariates, thus allowing for an observable trait in the graph, a similar strategy is adopted by roch2018generalized.
Though we emphasise the estimation of the stochastic blockmodel, there are several other methods available for estimation of the network formation model. These include the beta model of chatterjeed2011, in which the graph generation model consists of two model parameters $\nu_i, \nu_i$ possibly altered through some additional dyadic covariates $X^*_{ij}$
where $\tilde f$ is a link function. Alternatively one can consider the latent space model of Hoff2002LatentAnalysis which include latent positions on some unobserved manifold $\mathcal{M}^p$.
In each of these cases breza2023consistently illustrate consistent estimation rates in the $||{\hat \theta - \theta_0}||_{\infty} = \mathcal{O}_P\left(\sqrt{\frac{\log(n)}{n}}\right)$ with the use of aggregated relational data. Since this represents the coarsest datatype we expect similar rates to hold for subgraph sampling and respondent driven sampling. Though this rate is too slow for the to ignore the effect of the estimation of the graph model, in examples where one expect a high level of correlation among the outcomes it can be practical to use these methods.
Here we elaborate on the computation of a Z estimator. In general, an estimator may require specific implementation, we provide an illustrative example with logistic regression. Recall the characterization of the average estimating function $m_i(Y_i, \mathbf{a}, \mathbf{X}; \beta, \theta) = \mathbb{E}[\tilde m(Y_i, S_i(\mathbf{X}, G), V_i(\mathbf{a}, G); \beta)|\mathbf{Y}, \mathbf{a}, \mathbf{X}; \theta]$. Under this model, $P(Y_i = 1|S_i(\mathbf{X}, G), V_i(\mathbf{a}, G)) = \Lambda(\tilde h(S_i, V_i)^T \beta)$.
In order to compute the new estimating function, we need to be able to consider the distribution of the graph, conditional on the observed outcome $Y_i$. Specifically.
In a standard missing data problem, one would impute the missing covariates directly, however, due to the dependence through the graph, this can be very difficult to achieve in practice. However, it will be straightforward to sample from the graph model $P(G|\theta)$. Using a simple approach, we can compute the maximizer exploiting standard software methods using an EM algorithm dempster1977maximum, wu1983convergence. Suppose that we draw a sample of graphs from the generative model $\{G^{(l)}\}_{l = 1}^L \sim_{iid} P(G|\theta)$.
Let $w_i(Y_i, G; \beta)$ define the weight of an observation.
We next construct the EM algorithm as follows.
In practice, this allows for one to use standard solvers for the (M-step), after sampling a single time with the (E-step).
Additionally, one can include correlations across the observations $Y_i$ through the use of a generalized estimating equation approach. In other generalized linear models, additional assumptions may be required in order to model the full conditional distribution $P(Y_i|S_i(\mathbf{X}, G), V_i(\mathbf{a}, G); \beta)$ such as a dispersion component.
For many problems, the parameter of interest is a causal query conditional on the complete graph $G$ as described in section (ref). For example, one may care about the expected number of adoptions after seeding an individual in block $k$ v.s. block $k'$. In this section, we illustrate how to construct an estimate of the causal parameter $\Psi(\mathbf{a}| G)$ using our conditional model estimation procedure.
Let $\Psi(\mathbf{a}|\theta_0) = \mathbb{E}[ \Psi(\mathbf{a}|G)| \mathbf{a}, \mathbf{X}, \theta_0]$ be the average causal effect of policy $\mathbf{a}$ over all draws of the graph model $\theta_0$. We will establish conditions under which these two quantities are close to one another.
Recall the true conditional mean function $\mathbb{E}[Y|S_i = s, V_i = v] = h_0(s,v)$. Under a correctly specified conditional model, $h_0(s,v) = h(s,v; \beta_0)$, and $\Psi(\mathbf{a}|\theta_0) = \Psi(\mathbf{a}|\beta_0, \theta_0)$ where
In order to estimate $\Psi(\mathbf{a}| G)$ we plug-in the estimates for the mean model and network model $\Psi(\mathbf{a}|\hat \beta, \hat \theta)$. We next discuss the asymptotics of the plug-in estimate.
This lemma is essentially an application of the delta method, with the additional caveat that we estimate $\theta$ before the plug-in estimate. As before, this requires a fast estimate of the graph generative model parameter, but we add the slightly different assumption ((ref)) that the smoothness in the model class is over the conditional response models $\mathbb{E}[h(S_i, V_i; \beta) | \theta]$, rather than the estimating function $\tilde m(Y,S,V|\beta,\theta)$.
Convergence of the causal parameter to the average over graphs
As we have previously discussed, we can only hope to estimate $\Psi(\mathbf{a}|\theta_0)$ as we do not have access to the full graph $G$. We next introduce a simple conditions under which the parameter $\Psi(\mathbf{a}| G)$ is close to its average over draws of the graph $G \sim \theta_0$, $\Psi(\mathbf{a}| \theta_0)$.
The proof is a one-line application of McDiarmid's inequality. Previous related work such as breza2023consistently typically assume that such a quantity is consistent, however here we quantify the rate here. We next highlight an example;
Here we illustrate the optimal design approach for Z-estimators. In this example, the variance itself may depend on the a parameter $\beta$, and thus one can include a working candidate for the parameter $\beta'$. In general, one could also propose a feasible range of the working parameters $\beta'$ and consider the worst case variance in that range.
As an extension of our variance minimizing procedure, we can incorporate the uncertainty in our estimates of the model parameters. For instance, consider the following parametric bootstrap approach for estimating the model parameters of the stochastic blockmodel when using ARD.
For example, consider a scenario where we utilize the stochastic blockmodel and we collect ARD. Denote $\hat \theta = (\{\hat Z_i\}_{i = 1}^n, \mathbf{\hat P})$ the initial estimate of the model as computed from Lemma (ref). We can construct a sampling distribution of $\hat \theta^{(b)}$ using the following procedure. Let $X^*_{it}$ denote the ARD responses of the number of connections individual $i$ has to someone of trait $t$ and let $T_i \in \{0,1\}^T$ denote the trait memberships of the corresponding individuals.
This approach can work for any procedures which can allow for a sampling distribution of the model parameters $\{\hat \theta^{(b)}\}_{b = 1}^B$. For example baraff2016estimating considers a nonparametric bootstrap for respondent driven sampling.
In all such cases we would like to include thee uncertainty in $\hat \theta$ to the saturation assignment, we apply Algorithm (ref) (or Algorithm (ref)) to each of the $b$ draws. Using the distribution of variances obtained over the $b$ draws, one can compute average or upper confidence bounds on the variance. For example in the simulation illustrated in section (ref) we select saturations based on the 2-standard deviation upper confidence bound of the average variance across $b \in \{1,2,\dots, B\}$.
In this section, we introduce the proofs for the results in the main paper as well as the additional theoretical results presented in section (ref).
We first include a useful lemma for bounding the approximation of the error of the graphon model.
We now proceed with a the proof of the lemma.
Here we provide additional details with respect to several aspects of our methodology. We also include further details on several of the details for the implementation of competing methods in section (ref).
In our simulation setup in section (ref) we can also compute confidence intervals based on the regression $Y_i = \beta^T \mathbb{E}[\tilde h(S_i, V_i)] + \epsilon_i$ where we apply the Eicker-Huber-White sandwich estimator of the variance. We then compute the corresponding plug-in estimator of the variance using the covariates observed and Lemma (ref). Since the covariates in the true regression model behave like averages over the graph, we expect Lemma (ref) to hold and therefore the difference between the GATE for any one draw of the graph, and the true GATE is very small. We see in Figure (ref) that the coverage tends to be larger than the nominal 95%, though in general, due to model misspecification of the true-graph, there can be additional uncertainty due to the misspecification of the graph model. However, we see in this simple example that the coverage performs well with an off-the-shelf implementation.
We next consider an example using a local diffusion process. We suppose that seed nodes are placed at time $0$ and that outcomes are measured at time $T = 1$, allowing for diffusion to only take place to the immediate neighbors with a fixed probability $q$. In this case, for non-seed nodes the probability of infection is related to the total number of treated neighbors through the following link function. Under this model let $V_i \in \{0,1\}$ denote the exposure as to whether one of their neighbors have received the treatment, i.e. $V_i = I(\sum_{j}G_{ij}a_j > 0)$. Then
In this experiment, a single individual is seeded in each network. Our goal is to identify the best individuals in each of the network to seed and rank them by the expected variance of the estimator. We compare this to random seeding of individuals in the network as well as seeding by only the highest degree nodes. We use the networks constructed by the union of all connections of banerjee2019gossip. We construct estimates of the stochastic blockmodel as the partial data example using $K = 3$ in each case. We construct the traits using ARD responses based on number of connections with the following traits outlined in the Appendix in section (ref). We also include an alternative where a beta-model chatterjeed2011 is used in place of the SBM for the degree seeding where further details on estimation are included in section (ref). We then draw samples of the graph using the parametric bootstrap to obtain a resampled distribution of ARD $\{\mathbf{X}^{*(b)}\}_{b= 1}^{B}$ for $B = 1000$. We identify the optimal treatment block for each parameter according to section (ref). We simulate $1000$ draws of the draws in the diffusion process for each true, and plot the associated bias and RMSE of the seeding strategies in Figure (ref) with a true diffusion parameter $q = 0.2$.
In the full data case, the optimal strategy would be to seed the highest degree node in each of the networks and measure whether each of their neighbors are infected at time $T = 1$. However, this poses a problem for the stochastic blockmodel as we are essentially picking an outlier to seed, which is different than a typical member of the block over draws of the process. This can be corrected for using a model which accounts for degree heterogeneity, in our case, the beta model. In our optimal seeding strategy, we find that the RMSE is lower in both the degree optimized strategy with the beta model, as well as the block optimized strategy with the SBM, than even the full data version with a completely randomized allocation, hence highlighting the role of the interplay of the model of the graph and the experimental design. This behavior is observed in Figure (ref).
In this example, we consider a problem of optimal treatment assignment after the outcome model is estimated. We consider an example where an outcome model is estimated and transported to a new population. In this example we suppose that there is some benefit $\beta_1 > 0$ to receiving a treatment, and some smaller benefit based on the fraction of the neighbors treated $0 < \beta_2 < \beta_1$. We wish to assign treatments in a way that will maximize the expected outcome $\Psi(\mathbf{a}|G)$ for each network.
Where $q_i := \frac{1}{d_i}\sum_{j = 1}^n G_{ij}a_j$ denotes the normalized number of treated neighbors. We simulate the data with $\beta_0 = 1$, $\beta_1 = 1$ and $\beta_2 = 1/2$ with $\sigma_i \sim N(0,1)$. We choose this form of a response function since it will be simple to solve with an off the shelf mixed-integer programming approach using CVXR Fu2020CVXR:Optimization.
We suppose that in each example there is only a budget for $B \in \{10,20,40,80\}$ treatments for each of the villages. The goal is to maximize the overall expected outcome. We consider the following competing procedures. In this case, we suppose that we have a single pilot network where we can learn the model and the goal is to maximize the benefit on the remaining networks. We use the same gossip diffusion networks as in sections (ref) and (ref).
We compare the following seeding strategies.
Let $\mathbb{E}[Y_i|\mathbf{a}] = \beta_0 + \beta_1a_i + \beta_2(1 - a_i) \sum_{k' = 1}^K \hat P_{\hat k_i k'} n_{t, k}$ and let $n_{t, k} = \sum_{j : k_j = k} a_j$. Therefore, the objective function.
where $\zeta = \frac{1}{d_i}\sum_{i = 1}^n \mathbf{P}_{k_i, \cdot}$ and $\mathbf{n}_t = (n_{t,1}, n_{t,2}, \dots, n_{t, K})$. In general, given a conditional model, one may fine tune the optimization approach to the particular challenges of evaluating the optimal treatment allocation. We partition each network into $6$ blocks.
We plot the expected average outcome under each of the treatment allocations for the remaining $68$ networks after learning a model from the first pilot network. We repeat this for the total number of treatments $B \in \{10,20,40,80\}$.
In Figure (ref) we find that based on our method, we can achieve higher average outcomes than simple models based on the block positioning alone, emphasizing the importance of considering the potential outcome model when optimal targeting.
We can also replicate the results of beaman2021can's study on the evidence of pitplanting. They consider 3 measures of information diffusion. Firstly, if an individual has heard of pitplanting, second, if they know how to pitplant, and thirdly whether they adopt pitplanting in their practice. In order to control for one's position in the network, the authors consider the distance between the optimal seeds using two other targeting methods, simple diffusion, and geo-targeting as well as complex contagion. They then compare the increased odds of con
Again, we generate synthetic covariates and apply a stochastic blockmodel in order to estimate $K = 8$ blocks within each of the networks. We plot the coefficients for the connection to exactly $1$ seed, $2$ seeds and within radius $2$ of at least $1$ seed in Figure (ref). We note that we run the same regression as in beaman2021can, however, some since the full network data includes some additional noise top preserve anonymity, we do not have the exact same estimates of the coefficients as in their paper, however, the conclusions are substantively the same.
To aid in reputability, we include additional details regarding the implementation of our methods as well as competing methods.
Another common model utilized for random graph formation is the beta model coined by chatterjee2011random. Namely these are a class of models that can be learned based on their degree sequence. We consider a version where each node has an affinity parameter $\nu_i$ and the probability of connection between each pair of nodes is $P(G_{ij} = 1) = \nu_i\nu_j$. Let $\nu_n = \sum_{i = 1}^n \nu_i$ Therefore, $\mathbb{E}[d_i = d] = \sum_{j \not = i} P(G_{ij} = 1) = \nu_i(\nu_n - \nu_i)$. The set of parameters $\{\nu_i\}_{i = 1}^n$ can be estimated using an iterative solution to the fixed point equation: $$ \nu_i^{(t + 1)} = d_i/(\nu^{(t)}_n - \nu^{(t)}_i)$$
We utilize the measured traits to construct responses for ARD questions for each individual for the networks in banerjee2019gossip. The constructed ARD include traits which ask "How many people do you know ..."
For the estimation of the GATE using banerjee2013diffusion, we use Leiden clustering and denote the clusters traits. When replicating the results of beaman2021can, only a subset of nodes have available covariate. As was done in our examples with banerjee2013diffusion, we construct synthetic traits using the clusters observed from Leiden clustering for $K = 10$. ARD is then constructed based on the connections to nodes of each trait.
The two estimators we compare for estimation of the global average treatment effect are the difference in means estimator $\widehat{\tau}_{DM}$ and the Horvitz-Thompson estimator $\widehat{\tau}_{HT}$. Let $E_{i0}$ and $E_{i1}$ denote the events that all neighbours of $i$ are untreated (including $i$ themselves) and treated respectively.
In general, the Horvitz-Thompson estimator will be unbiased, however, it can often suffer from high variance for two reasons. Firstly, the probabilities of the events that all nodes are treated may be exceedingly low, inflating this variance, and also, relatively few nodes receive the exposures under which all of their neighbours are treated or none of them are.
In the case where the spillover effects are relatively mild, often a difference in means approach to the estimator is preferred. The effect of cluster randomization on the MSE of this estimator has been further studied in the complete network brennan2022cluster, viviano2020experimental.