The New England Journal of Statistics in Data Science

Network data encode not only information about individuals but also the connections among them. Such data exist in many settings: students who are friends with one another [4, 33]; intersecting streets [27]; interconnected brain regions [30]; and social network users who “follow” or “friend” one another [47]. Network data can be represented by a graph where nodes represent individuals and an edge between two nodes indicates that there is a connection between them. For example, Figure 1 visualizes a Caltech Facebook network, which was collected in September 2005 [42, 43]. The version of the data used in this paper contains information about 770 users and was retrieved from the Network Data Repository [38]. In Figure 1, the circles (nodes) represent Caltech Facebook users and a line between two nodes (an edge) signifies a Facebook friendship between two users.

The dependency structure in network data poses distinctive challenges and opportunities and thus motivates specialized theory and methodology. In this paper, we focus on network experimentation, a topic that has gained recent interest given the widespread experimentation happening in social network companies like LinkedIn, Facebook, and Twitter [17]. These businesses regularly develop updates and new features to their platforms. To decide whether a new feature or update should be deployed for all users, these companies conduct controlled experiments on their user networks to test a new feature and evaluate its causal impact [34, 47].

Current literature on network experimentation generally focuses on the setting where the experiment is conducted on a given network of a fixed number of individuals. Connections among individuals on the network are further assumed to be known and unchanged over the duration of the experiment. In addition, the experiment usually involves only two experimental conditions: the proposed new condition (i.e., the treatment) and the existing or baseline condition (i.e., the control). This type of experiment is also known as an A/B test, since two conditions A (control) and B (treatment) are compared to determine which one performs better in terms of some business metric. For example, the outcome of interest may be a measure of engagement like time spent on the platform, and the goal of the experiment may be to determine which condition maximizes this metric. In general, the experimental outcome is measured on experimental units, who are assigned to the experimental conditions. In the social network setting, the experimental units are often the platform users.

There are two stages in running an experiment: designing the experiment and analyzing the results. In the design stage, the experimenter determines the treatment assignment for the experimental units. Then in the analysis stage, data collected from the experiment are used to estimate the effect of the treatment on the outcome of interest. A popular estimand in the literature is the global treatment effect (sometimes alternatively referred to as the average treatment effect), which is the average difference in experimental outcomes when everyone in the population is assigned to the treatment versus when everyone is assigned to the control [9, 12]. A significantly positive global treatment effect (in the case where higher outcome values are desired) would motivate a business decision to deploy the treatment instead of the control to all users and vice versa. As such, an accurate estimate of the global treatment effect is an important objective to consider when developing experimental design and analysis methodologies. Another relevant concern in network experimentation is the estimation of the network effect: the influence of the network on the outcome of interest.

Classical experimental design and analysis methodologies treat individuals (or experimental units) as independent entities, in the sense that the potential outcome of each unit is assumed to be independent of the experimental conditions the other units are assigned. This is known as the Stable Unit Treatment Value Assumption (SUTVA) [11]. However, when experimental units are connected on a network, SUTVA is often violated. For example, consider an experiment where the treatment is designed to increase the time a user spends on the social network. Suppose Connie is assigned to control and her friend, Tracy, is assigned to treatment. Due to the treatment, suppose Tracy spends more time on the network. Connie may observe Tracy’s increased activity and interact with Tracy more, resulting in an increased amount of time spent on the platform despite being in the control group. Thus, the fact that Tracy is assigned to the treatment affects the potential outcome of Connie. As such, SUTVA no longer holds and the potential outcome of an individual depends not only on their own treatment assignment but also the treatment assignment of other individuals in the network. Consequently, in a network setting, the global treatment effect depends not only on the experimental conditions but also on the structure of the network. Therefore, specialized methods are required for the design and analysis of network experiments.

There are two major classes of approaches to the design and analysis of experiments on networks. The first base their methods on the exposure framework [2, 12, 16], in which experimental units (nodes) are classified into exposure groups where group members have the same level of exposure to the treatment. In a network experiment, even when two units are assigned to the same experimental condition, they may be exposed to different levels of the treatment (e.g., based on the treatment assignment of those around them) and thus belong to different exposure groups. For example, units i and j may belong to different exposure groups if all neighbors of unit i are assigned to treatment while all neighbors of unit j are assigned to control even when units i and j have the same treatment status. The literature often defines exposure groups based on local treatment assignment patterns of neighbors [2, 12, 39]. Based on the defined exposures, various causal effects, such as the global treatment effect and others [20, 41], can be defined as contrasts of the exposure groups’ potential outcomes [2]. Accordingly, estimators can be derived using inverse probability weighting [2, 12]. As an example, Gui et al. [16] define the treatment exposure as containing units such that (i) they are assigned to the treatment and (ii) the majority of their neighbors are also assigned to the treatment. Similarly, units are said to belong to the control exposure if (i) they are assigned to the control and (ii) the majority of their neighbors are also assigned to the control. The potential outcome of the treatment exposure group serves as a proxy for the outcome observed if the whole network is assigned to treatment. Similarly, the control exposure group is assumed to behave as if the whole network is assigned to control. The global treatment effect can then be estimated by taking the (inverse probability weighted) difference between the average outcome of the treatment exposure group and that of the control exposure group.

To achieve a precise estimator of the global treatment effect, each exposure group needs to contain a large number of units. In other words, the exposure-based analysis strategy requires that experimental units are mostly surrounded by neighbors assigned to the same experimental conditions as themselves. This motivates a popular design strategy in the network experimentation literature based on graph-cluster randomization [8, 16, 23, 44]. Graph-cluster randomization designs first partition the network into clusters, i.e., subgraphs of nodes that are densely connected within each cluster and sparsely connected between clusters. These clusters are then randomly assigned to treatment or control and all units in each cluster are assigned to the same experimental condition. This ensures that most units share the same experimental condition as their neighbors. Eckles et al. [12] illustrate via simulations that graph-cluster randomization, together with inverse-probability-weighted estimators, reduces bias in estimating the global treatment effect in experiments on networks.

The second class of approaches is based on assuming a statistical model for the experimental outcomes. There is a rich literature on network regression models for inference and prediction problems [1, 7, 22, 29, 31]. However, for network experiments, models need to characterize and quantify the influence of the treatment assignment to the outcome of interest, compared to the control, while taking into account the network influence. Parker et al. [35] and Koutra et al. [26] posit that an experimental condition can affect the outcome of a unit in two ways: (i) through the unit’s own treatment assignment and/or (ii) through the treatment assignment of the unit’s neighbors. These effects are assumed to be fixed and additive and then modeled via an ordinary linear regression model. Basse and Airoldi [5] and Pokhilko et al. [36] instead suggest that the experimental conditions only assert influence via the nodes’ own treatment assignment. Network interference is then modeled via the correlation of random errors based on the network structure, in a similar fashion to the conditional autoregressive (CAR) model from the spatial statistics literature [3, 10].

In the model-based framework, once the model is formalized, experimenters can select design criteria to reduce, for example, the variance of the model parameters and/or functions of them. Designs that optimize the posed design criteria can be found using exhaustive search [26, 35], random search [35] or other optimal design search methods [4, 36]. In the analysis stage, model parameters will be estimated by fitting the model to the observed data using methods such as least squares or maximum likelihood.

Of the two classes of approaches, the exposure framework seems more popular in the network experimentation literature thanks to its focus on the global treatment effect. When the main goal of the experiment is to accurately estimate the global treatment effect, the exposure framework only requires the experimenters to define the treatment and control exposures. For instance, in the Gui et al. [16] method the experimenter must determine what percentage of a unit’s neighbors must be assigned to treatment for it to be classified into the treatment exposure group. After that, graph-cluster randomization can be applied as a general design strategy to reduce the bias of the global treatment effect estimate. Nevertheless, there exist pitfalls associated with this class of approaches. Chin [9] points out that social networks are usually locally dense, making it hard for algorithms to partition the network into separate clusters. In addition, experimental outcomes on units not classified into the treatment or the control exposures will be wastefully discarded in the analysis stage [16], resulting in a low effective sample size [39].

On the contrary, if the generating model for the experimental outcome is correctly specified, the model-based framework will have several advantages over the exposure framework. First, data from all units participating in the experiment can be utilized when fitting the model in the analysis stage. Second, additional inference beyond the global treatment effect can be achieved via model parameters and functions of them. Finally, design selection can be calibrated according to the experimenters’ specific interests, which may sometimes be something other than the global treatment effect.

In this paper, we attempt to unify many of the model-based approaches by introducing the general additive network effect (GANE) model, in which the experimental outcome of a unit is modeled as an additive function of the effect of its own treatment assignment and the effects that come from other treated or controlled nodes in the network. The model not only encompasses (as special cases) many network experiment models from the literature, but it also flexibly extends the manner in which network effects are modeled. In particular, the network influence from the treatment and control groups can be modeled differently in terms of both size and functional form. Moreover, the model is specified such that its parameters are interpretable. We show that (quasi) maximum likelihood estimation is possible for a family of model specifications. The resulting (quasi) maximum likelihood estimators are then proven to be consistent and asymptotically normal, which facilitates familiar likelihood-based inference.

Existing work in the model-based direction concentrates on inference for the model parameters themselves, however, in many cases, the global treatment effect is the main quantity of interest. To expand the utility of the model-based framework, we illustrate how quantities such as the global treatment effect can be expressed as functions of the GANE model parameters, which permits inference via the delta method [45]. Similarly, via the model parameters, the experimenters can also build design criteria and select experimental designs that optimize these criteria. Overall, through the GANE model, we provide a generalized model-based framework for the design and analysis of experiments on networks.

We further propose a particular specification of the GANE model, the power-degree (POW-DEG) specification, where the network effects are modeled as a nonlinear function of the number of neighbors assigned to treatment or control. Using the power as an additional parameter, the POW-DEG model specification gains additional flexibility in capturing the network effect compared to the linear counterpart that has been proposed in the literature [35]. We investigate inferential properties of the POW-DEG specification as well as other existing specifications in the literature via simulations on real-life networks. We find that the POW-DEG specification yields good estimates of the true global treatment effect, even under model misspecification. By examining designs in the context of the POW-DEG specification, we also find that cluster randomization and balanced treatment assignment are not necessarily optimal for this model specification and other design strategies should therefore be explored in the future.

As with any model-based approach, the GANE model framework faces the challenge of model selection. When the GANE model is estimated via maximum likelihood, experimenters can use popular model selection criteria such as AIC, BIC, etc. to choose the model specification that best fits the data. However, in the design phase, before data have been collected, experimenters must use their domain knowledge to select a suitable model specification, with which they can build design criteria and select a design. In the absence of information concerning which model specification is appropriate, we suggest the use of the POW-DEG specification because of its flexibility and good performance in terms of global treatment effect estimation and inference. Nevertheless, model selection for both the design and analysis of experiments on networks remains an open area for future study.

The rest of the paper is structured as follows. In Section 2, we introduce the general additive network effect (GANE) model framework. We present the theoretical results for maximum likelihood estimation in Section 3. In Section 4 we report the results of several simulation studies used to investigate the design and inferential properties of a few particular specifications of the GANE model. In Section 5 we conclude with a summary and discussion of several extensions to this work.

Suppose that the experiment is on a network of n experimental units. Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ denote the network where $\mathcal{V}$ denotes the set of experimental units (i.e., nodes) and $\mathcal{E}$ denotes the set of connections among experimental units (i.e., edges). As in most of the network experiment literature, we assume that $\mathcal{G}$ is known and fixed throughout the experiment. Furthermore, $\mathcal{G}$ is assumed to be undirected and simple, and so $\mathcal{G}$ can be represented by an $n\times n$ adjacency matrix A where for $i,j=1,2,\dots ,n$ If ${A_{ij}}=1$, unit i and j are said to be neighbors. The number of neighbors of unit i, i.e., the sum of the ith row (or ith column) in the adjacency matrix, is called the degree of unit i, which we denote by ${k_{i}}$. We assume the treatment vs. control setting and let ${Z_{i}}$ be the binary indicator denoting the experimental condition of unit i, where ${Z_{i}}=1$ if unit i is assigned to treatment and ${Z_{i}}=0$ if unit i is assigned to control. Let ${Y_{i}}$ denote the experimental outcome of unit i and further let $\mathbf{D}=(\mathbf{A},\mathbf{Z},\mathbf{Y},\mathbf{X})$ denote the experimental data where Z and Y are the treatment assignment and outcome vectors respectively, A is the adjacency matrix defined above, and X is a possible $n\times p$ matrix containing the units’ covariates.

Parker et al. [35] introduce the linear network effect (LNE) model where ${\epsilon _{i}}\sim \mathcal{N}(0,1)$ independently for all $i=1,2,\dots ,n$. The model assumes that when a unit is assigned to the treatment group, the unit itself will experience an effect of size τ while exerting an effect of size ${\gamma _{T}}$ on each of its neighbors. If a unit is assigned to the control, it will exert an effect of size ${\gamma _{C}}$ on each of its neighbors. The outcome of a unit i is then the sum of the baseline μ, the effect of the treatment assignment, and the sum of the effects from its neighbors plus some random error ${\epsilon _{i}}$. In other words, the experimental outcome of unit i is a linear combination of which experimental condition it belongs to and the numbers of its neighbors that are assigned to treatment (treatment degree) and control (control degree). The model offers a straightforward parameter interpretation: μ is the expected baseline outcome, τ is the effect of the treatment assignment, and ${\gamma _{T}}$ and ${\gamma _{C}}$ quantify the effect of the neighbors’ treatment assignments. The model is linear in the parameters and can therefore be fit using ordinary least squares.

However, on a network, the outcome of a unit may not only be influenced by the number of neighbors it has. It may also be influenced by the outcomes of its neighbors [1, 40], or by neighbors of neighbors, etc. To enable experimenters to model the experimental outcomes more flexibly, we propose a more general model for the design and analysis of experiments on networks, which we call the general additive network effect (GANE) model: where ${f_{T,i}}$ models the effect of treated units on unit i and ${f_{C,i}}$ models the effect of controlled units. Hence, ${f_{T}}$ and ${f_{C}}$ are functions of the network (via the adjacency matrix A) and the treatment assignment vector Z. They may also be functions of the outcome vector Y, a possible covariate matrix X, and a parameter vector $\boldsymbol{\eta }$. By separating the network effect based on the source (treatment or control), the GANE model not only allows the network effect to come from both treated and controlled neighbors, but the network effect can also be of different forms and sizes, i.e., ${f_{T}}$ and ${f_{C}}$ may have different forms and/or may depend on different parameters (which are all contained in the vector $\boldsymbol{\eta }$). The functional forms of ${f_{T}}$ and ${f_{C}}$ may be determined by the experimenters’ domain knowledge or other considerations.

The GANE model has a structure similar to that of the LNE model (2.1) by Parker et al. [35], in which μ parametrizes the baseline outcome, τ parametrizes the effect of the treatment assignment to a unit and ${f_{T}}$ and ${f_{C}}$ respectively model the influence the unit receives from other treated and controlled units. Therefore, it retains the clear interpretability of Model (2.1) while increasing the flexibility with which the influence of other nodes is modeled.

The GANE model encompasses several existing models in the literature, a few of which are identified below.

The Linear Network Effect (LNE) Model: Model (2.1) is a GANE model with ${f_{T}}$ modeling the treatment degree and ${f_{C}}$ modeling the control degree: where $\boldsymbol{\eta }={({\gamma _{T}},{\gamma _{C}})^{\top }}$.

The Local Aggregate (LAG) Model: Advani and Malde [1] propose that the outcome of a unit may be influenced by the sum of its neighbors’ outcomes in cases such as (i) a person’s criminal behavior might depend on the number of crimes committed by their friends; (ii) the number of items purchased by a customer may depend on the number of items purchased by their friends; and (iii) the number of hours a student spends studying may depend on the number of hours their friends spend studying. Accordingly, under the GANE model, functions ${f_{T}}$ and ${f_{C}}$ can respectively model the sum of the treated and controlled neighbors’ outcomes: where $\boldsymbol{\eta }={({\rho _{T}},{\rho _{C}})^{\top }}$.

The Homophily (HOM) Model: In the social economics literature, a unit’s outcome is often modeled as a function of the average outcome of its neighbors [7, 31]. This is often referred to as the homophily effect [40] which is summarized by the conjecture that “you are the average of the people around you”, and a result of people’s desire to conform to social norms [1]. Gui et al. [16] formalize a model with homophily effect for the analysis of a network A/B test as follows: That is, besides the homophily effect, each treated node will exert an effect of size γ on each of its neighbors. This model can be reparameterized and written in the GANE model framework with In this case, $\boldsymbol{\eta }={({\rho _{T}},{\rho _{C}},{\gamma _{T}})^{\top }}$ with constraints ${\rho _{T}}={\rho _{C}}=\rho $ and ${\gamma _{T}}=\gamma $.

The above are existing models that have been proposed and discussed in the literature. We next propose a new specification under the GANE framework.

The Power-Degree (POW-DEG) Model: Let us continue with the example in Section 1.1 where the outcome of interest is the amount of time that users spend on a social network platform. Recall, since Tracy is assigned to treatment, Connie (in the control group) may spend more time on the platform simply because she is Tracy’s friend. However, Connie’s time increase due to her first treated friend may not necessarily be equal to the time increase due to her 100th treated friend. In particular, as the number of treated friends grows, the effect of an additional treated friend is likely to decrease. This intuition is similar to the law of diminishing marginal utility in economics [15].

When the number of neighbors is low, for example, in agricultural settings where units are plots placed on a lattice, it may be reasonable to model the network effects homogeneously as in the LNE specification (2.3). However, in the social network setting where a user can have hundreds to thousands of friends, this linearity assumption seems less likely. As a result, we suggest modeling the network effect as a non-linear function of the number of neighbors (i.e., the degree). In particular, we propose the following GANE model specification with where $\boldsymbol{\eta }={({\gamma _{T}},{\gamma _{C}},\lambda )^{\top }}$. We call this the power-degree (POW-DEG) specification because the network effects are modeled as powers of the treatment and control degrees. The power parameter λ serves to temper the growth of network effects as the treatment and control degrees increase, and so we expect that $0\lt \lambda \lt 1$. However, in the interest of ample flexibility, we do not make this assumption. We allow for the possibility that $\lambda \gt 1$ and also the possibility of the LNE specification (2.3) arising as a special case when $\lambda =1$.

Which estimation method is appropriate for fitting the GANE model depends on the specification of the functions ${f_{T}}$ and ${f_{C}}$. For instance, estimation and inference for the LNE specification (2.3) can be performed using ordinary least squares. However, when ${f_{T}}$ and ${f_{C}}$ are more complicated such as in specifications (2.4) or (2.5), the GANE model becomes (spatially) autoregressive, in which case maximum likelihood estimation is required.

As discussed in Section 1.2, compared to the exposure framework, one of the advantages of the model-based framework is that more quantities can be easily estimated, assuming they can be expressed as functions of the model parameters. Here, we review several such quantities that may be of interest to experimenters.

Global treatment effect (GTE): As noted already, an important quantity for business decision-making is the global treatment effect (GTE), which is defined as the difference in average outcomes when everyone in the network is assigned to the treatment versus when everyone is assigned to the control. Note that GTE measures the treatment effect at the global level, instead of the individual level, taking into account the structure of the network and possible network effects. In the GANE framework, this quantity can be expressed as where the subscripts on D indicate that the functions are evaluated when all the experimental units are assigned to treatment or control, respectively. Note that the expectations in the second equivalence are only necessary when ${f_{T}}$ and/or ${f_{C}}$ are functions of the outcome Y. As shown, the GTE can be expressed as a function of the model parameters. Moreover, GTE can be decomposed into two components, the direct and indirect treatment effects, which aid interpretation.

Direct treatment effect (DTE): Traditionally, the direct treatment effect with respect to a treatment assignment vector Z is defined as the average difference in expected outcomes when the treatment assignment of one unit changes while the others remain the same, i.e., where ${\mathbf{Z}_{-i}}$ denotes vector Z without the ith element. This definition is difficult to interpret and calculate, especially when the network effect functions are complicated. So instead, we define the direct treatment effect as the expected difference in outcomes when a node is assigned to the treatment versus when a node is assigned to control, keeping the network effects fixed. In the GANE model framework, the direct treatment effect is simply Our definition is clear, easy to interpret, easy to calculate, and does not depend on any specific treatment assignment vector. Hence, it can be used across all specifications of the GANE model.

Indirect treatment effect (ITE): Interest may also lie in quantifying the amount of the global treatment effect due to the network. The ITE is therefore defined as the difference between the global treatment effect and the direct treatment effect. In the GANE framework, this is Hence, the indirect treatment effect can also be interpreted as the difference between the network effect induced by the treatment versus that induced by the control.

Each of these quantities can be expressed as functions of the model parameters μ, τ, and $\boldsymbol{\eta }$. When the expectations of ${f_{T}}$ and ${f_{C}}$ are known, estimates and hypothesis tests for these quantities can be developed based on parametric inference associated with the model. With respect to hypothesis tests, the experimenters may be interested in one or more of the following hypotheses.

When the GANE model can be estimated using maximum likelihood estimation, since DTE, ITE, and GTE are functions of the model parameters, Hypotheses 1, 3, and 4 can be tested using Wald-type tests, and Hypothesis 2 can be tested using a likelihood ratio test.

To design an A/B test, the experimenter decides which units are assigned to the treatment group and which units are assigned to the control group. This corresponds to specifying the treatment assignment vector Z where each unit i is assigned ${Z_{i}}=0$ or 1. Interest lies in determining an optimal treatment assignment.

Under the GANE framework, experimenters can define design criteria related to the variance of the parameters’ estimators. For example, D-optimality [37] aims to minimize the determinant of the variance-covariance matrix and is therefore a popular choice to generally minimize the confidence region for the parameters. However, in the network experiment context, as discussed in Section 2.3, GTE (or possibly DTE or ITE) is the primary quantity of interest. Therefore, in the experimental design simulations in Section 4.4, we will focus on the variance of GTE estimates as our design criterion.

With the design criterion determined, experimenters can find good designs using methods such as exchange algorithms [26, 35], or random search [35]. Random search refers to the process in which the design is chosen by first randomly generating a large number of designs, and then choosing the design with the best-evaluated design criteria. Exchange algorithms [26, 35] instead take a greedy approach by iteratively changing the treatment assignment for each unit in the direction of optimizing the design criteria. Via simulation, Parker et al. [35] find that random search, despite its simplicity, yields nearly as good designs as the computationally less efficient exchange algorithm. We thus use random search as a design selection strategy in our simulations in Section 4.4.

We have discussed in Section 2.2 that different specifications of the GANE model may require different estimation techniques. In Section 2.3, we also mentioned that when maximum likelihood estimation is possible, different hypotheses can be tested using the maximum likelihood framework. Therefore, in this section, we discuss the maximum likelihood estimation for the GANE model.

In order to obtain the likelihood of the outcome vector Y, we consider the family of GANE specifications where the outcome ${Y_{i}}$ either (i) does not depend on neighboring outcomes, or (ii) depends linearly on neighboring outcomes. That is, where $\boldsymbol{\eta }={({\rho _{T}},{\rho _{C}},{\gamma _{T}},{\gamma _{C}},{\boldsymbol{\varphi }^{\top }})^{\top }}$ and ${W_{T,ij}}$ (or ${W_{C,ij}}$) is the ${(i,j)^{th}}$ element of pre-specified weight matrix ${\mathbf{W}_{T}}$ (or ${\mathbf{W}_{C}}$). The diagonals of these weight matrices are zero, i.e., ${W_{l,ii}}=0$ for $l\in \{T,C\}$ and $i=1,\dots ,n$. In addition, ${g_{T,i}}(\boldsymbol{\varphi })$ and ${g_{C,i}}(\boldsymbol{\varphi })$ are real-valued functions, possibly depending on the parameter vector $\boldsymbol{\varphi }$, the experimental data D, but not the outcome vector Y. We can see that Model (3.1) generalizes all model specifications discussed in Section 2.2, in which the outcome of an experiment may depend linearly on other unit’s outcomes and/or possibly nonlinearly on other covariates. Model (3.1), however, excludes cases where the outcome of unit i is dependent on a nonlinear function of the outcome vector Y, which complicates the maximum likelihood theory.

To perform estimation, we consider the matrix form of Model (3.1) as follows where ${\mathbf{G}_{T}}(\boldsymbol{\varphi })$ denotes the $n\times 1$ vector of ${g_{T,i}}(\boldsymbol{\varphi })$ values and ${\mathbf{G}_{C}}(\boldsymbol{\varphi })$ denotes the $n\times 1$ vector of ${g_{C,i}}(\boldsymbol{\varphi })$ values. Let $\mathbf{M}(\boldsymbol{\varphi })=[{\mathbf{1}_{n}}\hspace{5.69054pt}\mathbf{Z}\hspace{5.69054pt}{\mathbf{G}_{T}}(\boldsymbol{\varphi })\hspace{5.69054pt}{\mathbf{G}_{C}}(\boldsymbol{\varphi })]$ be the model matrix, which depends on the parameter vector $\boldsymbol{\varphi }$. Further, let $\boldsymbol{\beta }={(\mu ,\tau ,{\gamma _{T}},{\gamma _{C}})^{\top }}$ and $\boldsymbol{\rho }={({\rho _{T}},{\rho _{C}})^{\top }}$. The model may be rewritten, isolating for Y on the left-hand side, as follows where $\mathbf{S}(\boldsymbol{\rho })={\mathbf{I}_{n}}-{\rho _{T}}{\mathbf{W}_{T}}-{\rho _{C}}{\mathbf{W}_{C}}$. The expression in (3.3) is possible if and only if $\mathbf{S}(\boldsymbol{\rho })$ is invertible. Lemma 1 gives sufficient conditions on $\boldsymbol{\rho }$ so that $\mathbf{S}(\boldsymbol{\rho })$ is nonsingular. Although the condition is based on any matrix norm, in practice, we can use the popular spectral norm [19] to derive the constraints. The proof of Lemma 1 is given in Appendix A.1.

With Y expressed as in Equation (3.3) and assuming $\boldsymbol{\epsilon }\sim \mathcal{N}({\mathbf{0}_{n}},{\sigma ^{2}}{\mathbf{I}_{N}})$, the log-likelihood function for Y is where $\boldsymbol{\theta }={({\boldsymbol{\rho }^{\top }},{\boldsymbol{\beta }^{\top }},{\boldsymbol{\varphi }^{\top }},{\sigma ^{2}})^{\top }}$ is the vector of all model parameters. If the normality assumption is not made, then (3.4) becomes the quasi log-likelihood [28] and the estimators $\hat{\boldsymbol{\theta }}$ that maximize (3.4) are called the quasi maximum likelihood estimators.

To find the maximum likelihood estimates, we take the first order derivatives with respect to $\boldsymbol{\beta }$ and ${\sigma ^{2}}$ and equate them to zero to obtain Note that these are the solution of an ordinary least squares that regresses the transformed outcome variable $\mathbf{S}(\boldsymbol{\rho })\mathbf{Y}$ on the covariate matrix $\mathbf{M}(\boldsymbol{\varphi })$ when $\boldsymbol{\rho }$ and $\boldsymbol{\varphi }$ are known. Plugging this into the log-likelihood (3.4), we can obtain the profile log-likelihood Then, we can find $\hat{\boldsymbol{\varphi }}$ and $\hat{\boldsymbol{\rho }}$ by maximizing ${\ell _{\text{P}}}$ using a grid search on their respective parameter space or using numerical algorithms such as the Nelder-Mead method [32]. The maximum likelihood estimates $\hat{\boldsymbol{\beta }}$ and ${\hat{\sigma }^{2}}$ can be obtained by plugging $\hat{\boldsymbol{\varphi }}$ and $\hat{\boldsymbol{\rho }}$ into Equations (3.5) and (3.6).

Note that in (3.5), it is implicitly required that $M{(\boldsymbol{\varphi })^{\top }}M(\boldsymbol{\varphi })$ is invertible, i.e., that the columns of the model matrix $M(\boldsymbol{\varphi })$ are linearly independent. Although unlikely, multicollinearity may exist. For instance, in the POW-DEG specification (2.6), when the graph is fully connected (i.e., every node is connected with one another), or when the treatment and/or control degrees are the same for all nodes, the model matrix will have linearly dependent columns. It is thus important in the design stage to choose a design that ensures the model matrix has full rank.

Here, we study the behavior of the maximum likelihood estimators as the network size increases to infinity. We use the subscript n to denote the data for a given network size n. Model (3.3) then becomes where ${\mathbf{S}_{n}}(\boldsymbol{\rho })={\mathbf{I}_{n}}-{\rho _{T}}{\mathbf{W}_{Tn}}-{\rho _{C}}{\mathbf{W}_{Cn}}$. Let ${\boldsymbol{\theta }_{0}}={({\boldsymbol{\rho }_{0}^{\top }},{\boldsymbol{\beta }_{0}^{\top }},{\boldsymbol{\varphi }_{0}^{\top }},{\sigma _{0}^{2}})^{\top }}$ be the true parameter values. The consistency and asymptotic normality properties of the maximum likelihood estimators ${\hat{\boldsymbol{\theta }}_{n}}$ are given in Theorem 1 below.

The proof of Theorem 1 is given in Appendix A.3, following the ideas of Lee [28], treating ${\mathbf{S}_{n}}(\boldsymbol{\rho })$ and ${M_{n}}(\boldsymbol{\varphi })$ as non-stochastic for any given $\boldsymbol{\rho }$ and $\boldsymbol{\varphi }$. The random errors ${\epsilon _{i,n}}$ are assumed to be independent and identically distributed with mean zero and variance ${\sigma _{0}^{2}}$. When ${\epsilon _{i,n}}$ follow a normal distribution, ${\hat{\boldsymbol{\theta }}_{n}}$ is the maximum likelihood estimator (instead of a quasi maximum likelihood estimator), and we have ${V_{n}}(\boldsymbol{\theta })={J_{n}}(\boldsymbol{\theta })$ and

With the asymptotic normality result, inference for the parameters can be performed accordingly. The inference for the causal quantities given in Section 2.3 can then be carried out via the Delta method [45]. In particular, the global treatment effect for Model (3.1) is calculated as Using the Delta method, the variance of the GTE can be written as where $\mathbf{t}=\frac{\partial \hspace{2.5pt}\text{GTE}(\boldsymbol{\theta })}{\partial {\boldsymbol{\theta }^{\top }}}$. As $\text{DTE}(\boldsymbol{\theta })=\tau $ and $\text{ITE}(\boldsymbol{\theta })=\text{GTE}(\boldsymbol{\theta })-\tau $, inference for DTE and ITE can be derived in a similar manner.

We first investigate the asymptotic properties of the maximum likelihood estimates derived in Section 3. As the theoretical results concern the case where the treatment assignment vector Z is known, in this simulation, we fix a particular design where half of the nodes are randomly assigned to treatment and the other half are assigned to control.

First, we investigate the results for the POW-DEG specification (2.6) by generating outcomes on the given network (either the Caltech or UMich Facebook network) with the following parameter settings: $\boldsymbol{\beta }={(0,1,0.5,0.1)^{\top }}$ and $\sigma =1$. We further vary the power λ within the set $\{0.5,0.75,1,1.25\}$, where $\lambda =1$ corresponds to the LNE specification (2.3). With each combination of parameters, 1,000 runs are conducted where the outcomes are generated and the maximum likelihood estimates are calculated accordingly.

The distribution of the parameter and GTE estimates for the POW-DEG specification (2.6) are plotted in Figure 2. We can see that the distributions of all estimates are reasonably bell-shaped and symmetric and centered around the true values (dashed vertical lines) as is expected given the asymptotic theory. While the distribution of $\hat{\tau }$ remains the same under different values of λ, the variances of the other estimators decrease when λ increases. This is because the ranges of values within ${\mathbf{G}_{T}}$ and ${\mathbf{G}_{C}}$ in the model matrix M increase as λ increases, which in turn decreases the variance of the parameter estimates.

The coverage of 95% asymptotic confidence intervals and variances of the parameter estimates are given in Figure 3, where left axes correspond to variances and right axes correspond to coverage. The blue lines depict the asymptotic variances derived from $\mathbf{J}({\boldsymbol{\theta }_{0}})$ and the red lines depict the sample variance of the 1,000 parameter estimates. The agreement between these lines suggests that the asymptotic variances may be used reliably for inference. With respect to coverage, the coverage rates for the 95% confidence intervals are plotted as grey bars on the right axes and the dotted lines serve as a reference at 0.95. We can see that the obtained confidence intervals have the correct coverage. To summarize, the simulation corroborates the asymptotic theory and indicates that maximum likelihood procedures work as expected for the POW-DEG specification (2.6). Simulation results for a different set of parameter values are included in Section S1 of the Supplementary Material. These results suggest that when the network effect is small and the network size is moderate, consistent estimation of ${\gamma _{T}}$, ${\gamma _{C}}$, and λ is more difficult. However, the battery of simulations was also run on the UMich Facebook network, whose size ($n=3,749$) is almost 5 times that of the Caltech network, and we find that estimation of all parameters, whether the network and treatment effects are large or small, agrees with the asymptotic theory. These results are also available in Section S1 of the Supplementary Material.

We conducted a similar simulation study on the HOM specification (2.5) with $\mu =0$, $\tau =1$, ${\gamma _{T}}=0.5$, ${\rho _{T}}={\rho _{C}}=0.1$ and ${\sigma ^{2}}=1$. The results for both the Caltech Facebook network and the UMich Facebook network are shown in Figure 4. As with the POW-DEG (2.6) estimates, and in agreement with the likelihood theory, the distributions of these parameter estimates are bell-shaped and centered at the true values. Moreover, since the UMich Facebook network is larger, the variation in the estimates decreases, as expected. Notice that the true values of GTE are different for the two networks, even though all parameters used are the same. This illustrates how the true value of GTE depends not only on the parameters but also on the structure of the graph. Variances and confidence interval coverage are also plotted in Figure 4. As we would expect, the asymptotic variances are suitable for inference and the asymptotic confidence intervals have acceptable coverage. To demonstrate the generality of these findings we present additional simulation results for another set of parameter values in Section S1 of the Supplementary Material. The theory developed in Section 3 and the simulations presented here (for multiple GANE specifications, parameter values, and networks) demonstrate the general utility of maximum likelihood inference with GANE models.

As discussed in Section 2.3, under the GANE framework, we can test hypotheses about the DTE, SUTVA, the ITE and the GTE. In particular, testing DTE = 0 is equivalent to testing ${H_{01}}:\tau =0$; testing whether SUTVA is satisfied is equivalent to testing ${H_{02}}:{f_{T}}={f_{C}}=0$; the null hypothesis for testing the indirect treatment effect is ${H_{03}}:\text{ITE}=0$; and the null hypothesis for testing the global treatment effect is ${H_{04}}:\text{GTE}=0$. In the maximum likelihood framework, Hypotheses 1, 3, and 4 can be tested using Wald-type tests and Hypothesis 2 can be tested with a likelihood ratio test.

We study the characteristics of these tests via simulation. Again, as the design Z is treated as fixed in our analysis, we randomly pick a design where half of the nodes are assigned to treatment and the other half are assigned to control. The parameters of the POW-DEG specification (2.6) are set at $\boldsymbol{\beta }={(0,1,0.5,0.1)^{\top }}$, $\sigma =1$ and $\lambda \in \{0.5,0.75,1.00,1.25\}$. Separate simulations are conducted to investigate each of the four hypothesis tests. For each simulation, values of certain parameters in $\boldsymbol{\beta }$ vary while the others stay as stated. In particular, in the simulation for Hypothesis 1, τ varies in the range $[0,1]$; in the simulation for Hypothesis 2, ${\gamma _{T}}={\gamma _{C}}$ and their values vary in the range $[0,0.05]$; for Hypothesis 3, ${\gamma _{C}}$ is fixed at 0.1 and ${\gamma _{T}}-{\gamma _{C}}$ varies in the range $[0,0.5]$. Hypothesis 4 with ${H_{04}}:\text{GTE}=0$ is also tested within each of the three simulations (with different λ) and the results are aggregated over different values of GTE, which corresponds to different parameter combinations. All tests are done at a 5% significance level and 1,000 runs are conducted for each parameter combination. The results are presented in Figure 5. The dotted horizontal line serves as a reference at the 5% level.

As expected, the rejection rates for each test increase as the respective parameter values depart from their null values. Moreover, tests for ${H_{01}}:\tau =0$ seem to behave similarly over different values of λ. This is consistent with the model estimation results in Figure 2 where the variances for $\hat{\tau }$ and $\widehat{\text{GTE}}$ look similar over different values of λ while the variances for $\hat{\gamma }$’s decrease as λ increases. We also remark that the results at null values deviate slightly from the nominal 5% level. This can be attributed to the use of asymptotic (inexact) variances in these tests. Although we do not include the results for the UMich Facebook network, given its size and given the results from Section 4.1, we expect similar results to those presented here for the Caltech Facebook network.

We conducted similar simulations for the HOM specification (2.5) with $\mu =0$, $\tau =1$, ${\gamma _{T}}=0.5$, ${\rho _{T}}={\rho _{C}}=0.1$, and varying τ, ${\gamma _{T}}={\rho _{C}}$ and ${\gamma _{T}}-{\rho _{C}}$ in different simulations for different hypothesis tests. The results are included in Section S2 of the Supplementary Material. It can be noted that the results are similar in both networks, except for different values of τ, signifying that this is an important parameter for the HOM specification (2.5). We also note that the rejection rates for Hypothesis 3 always stay at $100\% $ even at ${\gamma _{T}}={\rho _{C}}$ (i.e., when the scaling coefficients for ${f_{T}}$ and ${f_{C}}$ are equal). This shows that the indirect effect is not only affected by the sizes of the network effects, but also by the functional forms of ${f_{T}}$ and ${f_{C}}$. These are different in the HOM specification (2.5).

The simulations in Sections 4.1 and 4.2 explore properties of maximum likelihood inference for different GANE specifications when they are correctly specified. In this section, we further investigate properties of these specifications under model misspecification. The specifications considered here are (i) the SUTVA specification, in which network effects do not exist and ${f_{T}}={f_{C}}=0$; (ii) the linear network effect (LNE) specification in (2.3); (iii) the POW-DEG specification in (2.6); (iv) the local aggregate (LAG) specification in (2.4); and (v) the homophily (HOM) specification in (2.5).

In this simulation, on the Caltech Facebook network, outcomes are generated 1,000 times for each of the listed model specifications. The data are then fitted using each of the five model specifications. We use global treatment effect (GTE) estimation and its inference results to compare performance among specifications because the GTE is generally of primary interest. To make the comparison fair, parameters for each model specification are chosen such that the true global treatment effect (GTE) is fixed at 2.0 and the average outcome variance is 1.0 in all data-generating scenarios. The exact parameter values for each specification are provided in Table 2.

Results of the simulation are plotted as heatmaps in Figure 6. The columns correspond to outcome-generating models and the rows correspond to estimating models. The top left panel shows the log ratio of the average estimated GTE to the true GTE. The desired value is 0, which is colored white. Red represents overestimation and blue represents underestimation. We see that all specifications can estimate the SUTVA specification well because it is nested within all GANE specifications. The POW-DEG specification seems to provide estimates with the lowest bias, even under model misspecification.

The top right panel shows the variances of the GTE estimates where white represents low variances and dark green represents high variances. We can see that the SUTVA and HOM (2.5) specifications provide the lowest variances while the highest variances come from the LAG specification (2.4). Both the POW-DEG (2.6) and the LNE (2.3) specifications provide reasonably low variances.

The bottom left panel shows the coverage rate of 95% confidence intervals for the GTE constructed by each estimating model. The results show that LNE (2.3), POW-DEG (2.6) and LAG (2.4) specifications have high coverage rates while the HOM specification (2.5) has lower coverage rates and the SUTVA specification has the worst. This is because the SUTVA specification does not capture the network effects introduced by other specifications.

Finally, on the bottom right, the model selection results by AIC are presented. Green represents high selection rates while white represents low selection rates. AIC works well as it selects the correct model specification most of the time, which is shown by the green diagonal. This supports the use of likelihood-based model selection criteria such as AIC for the GANE framework. Furthermore, it can be seen that POW-DEG specification (2.6) is selected fairly often no matter the data-generating model, which suggests that it fits the data reasonably well even under model misspecification.

As illustrated in Figure 6, the POW-DEG specification (2.6) is the only one that performs well in each dimension. This illustrates the flexibility of the POW-DEG model to capture a variety of network effects. Hence, we advocate its use generally, especially when there is no prior information or preference for another specification.

Given the suggestion for the general use of the POW-DEG specification (2.6) in Section 4.3, in this section, we examine the characteristics of good designs in the context of this specification. Identifying such characteristics will provide useful insights to consider when designing the experiment. To do so, we randomly generate 10,000 designs for the Caltech Facebook network as follows: (i) take a random draw m from the discrete uniform distribution on $[1,n-1]$; (ii) randomly select m nodes and assign them to treatment, and assign the remaining $n-m$ nodes to control.

As GTE is the primary focus for many experiments on networks, in this simulation we use the variance of the estimated global treatment effect $\mathrm{Var}[\widehat{\text{GTE}}]$ as the design criterion to evaluate the designs. Good designs are ones that give lower values of the design criterion. The criterion is evaluated based on the Fisher information matrix of the POW-DEG specification (2.6) with parameters $\mu =0$, $\tau =1$, $\lambda \in \{0.50,0.75,1.00,1.25\}$. We also vary ${\gamma _{T}}$ and ${\gamma _{C}}$ in three settings: (i) ${\gamma _{T}}=0.5\gt {\gamma _{C}}=0.1$; (ii) ${\gamma _{T}}=0.5={\gamma _{C}}$; (iii) ${\gamma _{T}}=0.1\lt {\gamma _{C}}=0.5$.

For each parameter combination, the best and worst 500 (i.e., 5%) designs with respect to $\mathrm{Var}[\widehat{\text{GTE}}]$ are recorded. We investigate three characteristics of these designs using three metrics. First, as Bowers et al. [6] point out, in the presence of network interference, balanced allocation may not be optimal. As such, we examine the percentage of treated nodes in each of the designs. Second, Parker et al. [35] observe that for certain design criteria, it is better to assign nodes with high degrees to treatment and nodes with low degrees to control. Hence we calculate the difference in average degrees of treated and controlled nodes to verify this observation. Third, if the best designs are clustered as graph-cluster randomization algorithms suggest, then nodes will be more likely to be surrounded by neighbors who share the same treatment assignment as themselves. This means that for each node, the percentage of neighbors having the same treatment assignment will be high. However, with our two-step design-generating procedure, we are considering both balanced and unbalanced designs. Suppose a design has b% of treated nodes, then the nodes will be expected to have b% of neighbors treated and $(100-b)$% of neighbors controlled. Nodes in a clustered design will have a high percentage of similarly treated neighbors regardless of such expectation. Therefore, we use the percentage of treated neighbors in each design as a threshold and calculate the percentage of nodes in the design that have the proportion of neighbors sharing the same treatment with themselves higher than this threshold.

Figure 7 shows characteristics of the best 500, the worst 500, and 500 randomly selected designs from the 10,000 that were generated. First, let us look at the best 5% designs. We see that the designs are more balanced at higher values of λ. At smaller values of λ, when ${\gamma _{T}}$ is greater, we need more treated nodes and when ${\gamma _{C}}$ is lower, we need more controlled nodes. Thus, in order to decide which treatment allocation to use, we need to know which of ${\gamma _{T}}$ or ${\gamma _{C}}$ is greater. This information may or may not be available in practice. Second, the differences in average degrees between treated nodes and controlled nodes distribute evenly around zero, suggesting that there is no preference in terms of degree when assigning treatments to the nodes. Last, in the third panel, we only observe a slight deviance of the boxplots from the 0.5 reference line at low values of λ. This suggests that graph cluster randomization is not very effective, at least in this specific setting. In contrast, the random designs have all the characteristics randomly distributed around the reference lines, while the worst designs will assign almost all of the nodes to either treatment or control.

Table 3 shows the efficiency gained (lost) by selecting the best (worst) 500 designs out of 10,000 random designs relative to a randomly selected design. We can see that the efficiency gained by the best 500 designs increases when network effects are present, i.e., when $ITE\ne 0$ or ${\gamma _{T}}\ne {\gamma _{C}}$ for the POW-DEG specification (2.6). It also increases when the sizes of the network effects increase, and when λ increases. Hence, even though a randomly generated design may be suitable when there is no network effect, when some network effect is present, other design selection techniques are necessary. We perform an analogous investigation with the UMich Facebook network, the results of which are provided in Section S3 of the Supplementary Material. The findings are similar to those presented here for the Caltech Facebook network, but we also notice that with a larger network, the efficiency of this random search procedure decreased. This suggests that when there are many experimental units, a balanced randomized design may still be efficient.