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.
42,510 characters · 6 sections · 45 citation commands
Higher-Order Causal Message Passing for Experimentation with Complex Interference
Randomized experiments are widely recognized as a reliable method in data-driven decision-making for determining the causal effects of new interventions, such as medical treatments or upgrades of market products. The conventional approach involves administering the new treatment to a randomly selected subset of the observation units (e.g., patients, products, or geographical areas), referred to as the treatment group, and comparing their outcomes with those units who received no treatment, the control group. However, the validity of these methods substantially relies on the assumption that treating a group of units does not interfere with the outcomes of the control units, known as the Stable Unit Treatment Value Assumption (SUTVA) cox1958planning, rubin1978bayesian, manski1990nonparametric,imbens2015causal,sussman2017elements.
In many social science and online marketplace scenarios, treating one unit impacts not only its outcome but also the outcomes of units that directly or indirectly interact with the treated unit bond201261, blake2014why, holtz2020reducing, johari2022experimental,bright2022reducing. This interference of treatments and outcomes makes estimating the causal effect of the treatment particularly challenging. Considering the network of interactions, when a unit is treated, its interactions with neighboring units lead to subsequent changes in their outcomes. These interactions continue over the experimental time horizon and may display complex behaviors. For example, as the treatment is expanded to a larger population, the interference effect may intensify or diminish. This necessitates efficient data usage and robust estimators to capture and adapt to such intricacies.
Given the complexity of analyzing interference phenomena, research on network interference often relies on a series of simplifying assumptions. One common assumption is to ignore variations over time and assume outcomes are observed at equilibrium, which discards valuable information before the system reaches equilibrium. To reduce the complexity of the analysis, further assumptions are imposed on the nature and level of interference choi2017estimation,cortez2022staggered, li2022network, such as the neighborhood interference assumption or assumptions on the maximum degree of the network. Additionally, a frequently made assumption to help estimate treatment effects is that the interference network is observed chen2024optimized, agarwal2022network, jia2024clustered, which is impractical in some settings, such as under pervasive interference. For example, in large-scale online platforms, units may interact through competing platforms, making it difficult to account for all sources of interference. Our aim in this paper is to relax these assumptions.
The impact of network interference can be intricate, particularly when considering interactions among units over time. For example, applying the treatment to one unit can have spillover effects on some control units, or one unit’s outcome can directly exert peer effects on other units’ outcomes. Simultaneously, treatments with long-lasting effects can have carryover effects to future time periods, and units’ outcomes can be serially correlated or have autocorrelation over time. Consequently, whenever SUTVA fails to hold, the number of potential outcomes grows exponentially with the population size and the time horizon of the experiment. This renders the estimation of causal effects under general interference structures impossible due to non-identifiability challenges manski2013identification,aronow2017estimating,basse2018limitations, karwa2018systematic,forastiere2022estimating.
Recently, shirani2024causal introduced a new framework called Causal Message-Passing (CMP) to address the challenge of causal effect estimation under unobserved pervasive interference. Their methodology relies on observing outcomes over time and is rooted in statistical physics mezard1986spin,mezard2009information and approximate message passing (AMP) donoho2009message,bayati2011dynamics from high dimensional statistics. Instead of investigating the complex relationships among units, which requires knowledge of the network, CMP focuses on the dynamics of one-dimensional quantities, such as the sample mean and sample variance of units' outcomes over time. These one-dimensional equations, also known as state evolution equations, can help track how the administered intervention propagates through the network of units over time, which enables the estimation of counterfactual scenarios. However, it remains underexplored how to use state evolution to estimate causal effects.
In this work, we propose to utilize machine learning to learn a mapping that updates key parameters of the distribution of outcomes over time for causal effect estimation. This is achieved by introducing a set of non-linear feature functions that act on the observed outcomes, creating a “basis” for the learning task. By training a properly designed machine learning model on this extracted basis, we estimate the Total Treatment Effect (TTE), also known as the Global Treatment Effect (GTE) or Global Average Treatment Effect (GATE), which measures the causal effect of altering the treatment scenario from treating no one to treating everyone. The result is a family of estimators that allow one to extract more information from the experimental data, thereby ensuring efficient use of the data.
To be more specific, this work builds on the foundation established by shirani2024causal, extending their method in two directions by introducing Higher-Order Causal Message Passing (HO-CMP) algorithms. First, HO-CMP incorporates higher-order moments of unit outcomes, unlike shirani2024causal's approach, which only employs the first moments for estimation. Second, while shirani2024causal focus solely on two-stage experiments with two different probabilities of treatment, our work leverages the additional data provided by having more than two experimental stages with multiple probabilities of treatment. Thus, our work aligns with the common practice in the tech industry of rolling out treatments through a sequence of experiments kohavi2020trustworthy.
We then validate the performance of HO-CMP by simulating multiple experimental settings, encompassing both linear and non-linear outcome specifications and various types of interference, such as synthetic random geometric networks and real-world networks. Specifically, we introduce a Non-LinearInMeans outcome specification, where the spillover effect is non-monotone in the fraction of treated neighbors; as an example of a complex treatment effect structure, we demonstrate how HO-CMP successfully estimates the total treatment effect by effectively utilizing higher-order moments of unit outcomes.
Simulating the experiments also allows us to calculate the ground truth value of the TTE, which remains unknown in real experiments, enabling us to compare the performance of HO-CMP to the ground truth TTE. Additionally, we benchmark HO-CMP against standard approaches such as difference-in-means and Horvitz-Thompson estimators, a recent technique of cortez2022staggered, and a first-order CMP estimation, like the one by shirani2024causal. We emphasize that a large body of recent estimators, e.g., jia2024clustered, requires knowledge of the interference network and is not applicable in our setting. The results showcase HO-CMP outperforming the benchmarks in estimating the TTE over time and its flexibility to cover different outcome specifications and interference structures.
\paragraph{Related causal inference literature.} The primary objective of research on causal inference in the context of network interference is to estimate causal effects while relaxing SUTVA. For this purpose, various assumptions and methods have been proposed. We briefly discuss the predominant ones.
A common approach to relax SUTVA is partial interference. Under this assumption, units are divided into disjoint clusters and interference is assumed only within the same cluster sobel2006randomized,rosenbaum2007interference,hudgens2008toward,tchetgen2012causal,liu2014large,kang2016peer,viviano2020policy,bhattacharya2020causal,qu2021efficient,auerbach2021local,candogan2023correlated,ugander2023randomized. When interference extends across clusters, standard estimators become biased. To address this, eckles2016design propose a cluster-randomized approach that randomizes treatment assignment across clusters, reducing bias. However, it requires knowledge of the clusters.
The other assumption to replace SUTVA is the Neighborhood Interference Assumption (NIA). NIA states that outcomes are only influenced by the treatments of neighboring units in the network. This assumption is commonly imposed in the literature that relaxes the SUTVA sussman2017elements. Some recent studies combine the NIA with the availability of either a fully or partially observed interference structure leung2020treatment, viviano2020experimental, agarwal2022network, belloni2022neighborhood, li2022random. Without prior knowledge of the interference structure, cortez2022staggered consider low-degree polynomial interactions among units in the network. leung2022causal also introduces a weaker version of the NIA, where the interference between two units located far away from each other is allowed to be nonzero, but negligible.
Another approach is to facilitate the estimation of causal effects by setting restrictions on the network structure chin2018central,jagadeesan2020designs,wang2020design,li2022network,agarwal2022network,jagadeesan2020designs,leung2022causal. These restrictions include bounding the largest node degree of the interference graph, limiting the degree of the dependency graph, observing specific patterns in the network, locally constrained interference structures, and restricting the topology of the interference network.
Driven by applications in marketplace platforms and two-sided marketplaces, several recent works have examined specific interference patterns holtz2020reducing,wager2021experimenting,munro2021treatment,johari2022experimental,harshaw2022designEC,farias2022markovian,bright2022reducing,farias2023correcting. For example, farias2022markovian study experiments in Markovian systems where interference effects propagate through constraints like limited inventory.
From another perspective, most of the existing literature on network interference focuses on the case of single-time point observation hudgens2008toward,aronow2017estimating,basse2019randomization,jackson2020adjusting,savje2021average. These studies have provided insightful results on spatial interference effects, but they often overlook temporal variations of the treatment effect. Recently, there has been a shift to consider settings with multiple-time observations li2022network,boyarsky2023modeling. However, the problem of considering the dynamics of units' outcomes remains understudied arkhangelsky2023causal.
Consider a system of $N$ units indexed by $i \in [N]:= \{1, \cdots, N\}$ subject to a randomized experiment. The units are observed over a time horizon of $T+1$ periods and for each $t \in \{0, 1, \ldots, T\}$, we let $\treatment{i}{t}$ denote the treatment status of unit $i$ during time period $t$. For simplicity, we consider a Bernoulli randomized design such that $\treatment{i}{t} \sim \text{Bernoulli}(\pi_t)$. That is, at time $t$ unit $i$ receives the treatment with a probability of $\pi_t$, corresponding to $\treatment{i}{t} = 1$. Otherwise, unit $i$ belongs to the control group and $\treatment{i}{t} = 0$. In this context, we collectively define $\bm{\pi} = (\pi_0, \pi_1, \ldots, \pi_T)$ as the experimental design. Then, following the potential outcome framework imbens2015causal, let $\outcome{i}{t}{}(\Mtreatment{}{})$ represent the potential outcome of unit $i$ at time $t$, where $\Mtreatment{}{}$ denotes the entire treatment allocation matrix, with $\treatment{i}{t}$ as the entry in row $t$ and column $i$.
Administering the treatment of unit $i$ at time $t$ according to $\Otreatment{i}{t}$ (as one realization of the random variable $\treatment{i}{t}$), we use $\MOtreatment$ (as one realization of $\Mtreatment{}{}$) to show the matrix that captures the treatments of all units throughout the experiment; accordingly, we let $\Ooutcome{i}{t}{} = \outcome{i}{t}{}(\Mtreatment{}{}=\MOtreatment)$ be the observed outcome of unit $i$ at time $t$ under the treatment assignment $\MOtreatment$:
Observing $(\MOtreatment, \MOoutcome)$, we are interested in estimating the TTE of the intervention, defined as below:
where $\bm{1}$ and $\bm{0}$ are matrices of all $1$ and all $0$ of appropriate dimensions (in this case, $T+1$ by $N$). Intuitively, the TTE measures the average effect of changing the treatment for the entire population. This is a common estimand in the network interference literature and provides important insights into the efficacy of the treatment for decision-makers jia2024clustered,chen2024optimized,viviano2023causal,yu2022estimating,cortez2022staggered.
Deriving a practical and efficient estimator for the TTE is challenging due to the fact that we can observe the population only under one treatment scenario holland1986statistics. Indeed, in Eq. (ref), we can observe at most one of $\outcome{i}{t}{\bm{1}}$ or $\outcome{i}{t}{\bm{0}}$, and often, neither.\footnote{We can only observe $\MOoutcome$ for one of exponentially many realizations of $\MOtreatment$.} In the following sections, we address this challenge by proposing a new class of estimators grounded in the CMP framework. These estimators rely on the efficient use of experimental data, $\MOoutcome$ and $\MOtreatment$, yielding accurate causal estimation under unknown network interference.
In this section we provide a summary of the outcome specification and results of shirani2024causal that we utilize in the remaining. For $t = 0, 1, \ldots, T-1$, we let $\outcomeg{t}{}: \RR \times \RR^{T+1} \mapsto \RR$ be an unknown measurable function. We also use $\Vtreatment{i}{} = \big(\treatment{i}{0}, \ldots, \treatment{i}{T}\big)^\top$ to denote the treatment assignment of unit $i$ during the experiment. Accordingly, the treatment allocation matrix $\Mtreatment{}{}$ is a $T+1$ by $N$ matrix with columns equal to $\Vtreatment{i}{}$. Given potential outcomes $\outcome{j}{t}{}(\Mtreatment{}{})$ at time $t$ and $j\in [N]$, their outcomes in time period $t+1$ are specified by
where $\IMatGl{ij}{}$ quantifies the impact of unit $j$ on unit $i$ at time $t$ and $\noise{i}{t}$ is a zero-mean Gaussian noise with a variance of $\sigma_e^2$, accounting for measurement errors. In addition, we let $\IMatG{} = [\IMatGl{ij}{}]_{i,j\in [N]}$ and refer to it as the interference matrix. Then, according to Eq. (ref), the function $\outcomeg{t}{}$ captures the impact of past outcomes and treatment assignments of other units on the current outcome of unit $i$.
Now, fixing $t$, we define
Then, as shown by shirani2024causal, whenever the elements of the interference matrix $\IMatGl{ij}{}$ are i.i.d. Gaussian random variables with mean $\mu/N$ and variance $\sigma^2/N$, under mild moment conditions on initial values $\outcome{i}{0}{}$, we have
where $Z_t \sim \cN(0,1)$ is independent from $\Vtreatment{}{} \sim \text{Bernoulli}(\bm \pi)$ (that is, $\treatment{}{t} \sim \text{Bernoulli}(\pi_t)$ and $\Vtreatment{}{} = (\treatment{}{0}, \treatment{}{1}, \ldots, \treatment{}{T})^\top$) and the equalities hold almost surely. We note that the theory behind this result is rooted in the AMP literature, going back to bolthausen2014iterative,bayati2011dynamics. However, as shirani2024causal note, there is a major distinction between the AMP literature and the above setting: in the AMP literature, the matrix $\IMatG{}$ is observed, and the aim is to construct proper functions $\outcomeg{t}{}$ for a completely different objective, which is studying the high-dimensional asymptotics of first-order algorithms. However, in the current context, the matrix $\IMatG{}$ and functions $\outcomeg{t}{}$ are unknown and the goal is to estimate them.
Considering Eq. (ref), the equations in (ref) determine the dynamics of the sample mean and sample variance of unit outcomes over time in large sample asymptotics, and are denoted by the State Evolution (SE) equations of the experiment shirani2024causal. In the next section, we present an efficient algorithm to learn the state evolution dynamics outlined in Eq. (ref). This method enables us to accurately estimate the TTE defined in Eq. (ref) and its corresponding confidence interval.
In this section, we introduce Higher-order Causal Message-passing (HO-CMP) for estimating the TTE over the entire time horizon of the experiment. Briefly speaking, HO-CMP directly estimates the update function in the state evolution equations (ref), thereby estimating counterfactual quantities while accounting for the impact of unknown network interference. To this end, by Eqs. (ref) and (ref), we rewrite the TTE as the difference of the sample means in the large limits: $$ \TTE_t = \nu_t(\bm{1}) - \nu_t(\bm{0}). $$ That means the problem of estimating the TTE is equivalent to estimating $\nu_t(\bm{1})$ and $\nu_t(\bm{0})$ using the observed data, denoted by $(\MOtreatment, \MOoutcome)$. On the other hand, considering the state evolution equations in (ref), the system state at time $t+1$, denoted by $(\nu_{t+1}(\MOtreatment), \rho_{t+1}(\MOtreatment)^2)$, is a (nonlinear) function of the system state distribution at time $t$, characterized by $(\nu_{t}(\MOtreatment), \rho_{t}(\MOtreatment)^2)$ and $\Vtreatment{}{}$, encompassing the sample mean and variance of observed outcomes as well as the design of the experiment. However, because the exact functional form and parameters of equations in (ref) are unknown, one cannot directly apply the SE to track the evolution of states. Therefore, we propose to estimate the unknown update functions in SE equations, utilizing the observed data $(\MOtreatment, \MOoutcome)$. For this purpose, we fix the treatment assignment matrix $\MOtreatment$ and define
In addition, let $\Vec{\phi} = (\phi_k)_{k\in[K]}$ be a prespecified vector of measurable feature functions of current estimates of the sample mean $\hat{\nu}_t(\MOtreatment)$, sample variance $\hat{\rho}_t(\MOtreatment)^2$, and the design $\MOtreatment$. We define $\bm{x}_t$ to represent the feature vector as follows:
Then, we formally propose learning the mapping $f_{\bm{\theta}}(\cdot)$ defined by,
We summarize the method in Algorithm (ref). Note in our experiment design we begin with all units under control by setting $\pi_0 = 0$, meaning no units receive treatment in period $0$. Additionally, to avoid non-identifiability issues, the experiment requires at least two stages, which corresponds to having at least two distinct values in the set $\{\pi_1, \ldots, \pi_T\}$.
The proposed HO-CMP method encompasses a rich family of estimators, offering flexibility through the selection of feature functions $\{\phi_k\}_{k \in [K]}$ and model $f_{\bm{\theta}}(\cdot)$. Specifically, incorporating proper feature (basis) functions, with examples shown in Table (ref), facilitates the extraction of informative patterns for learning the unknown nonlinear dynamics of the system throughout the experiment. In practice, one could choose these basis functions based on heuristics, domain knowledge, and prior information about the dynamics.
Specifically, in this paper, we consider the following estimators, as summarized in Table (ref).
FO-CMP (First-Order Causal Message-Passing): This corresponds to the simple setting where $\nu_{t+1}(\MOtreatment)$ is assumed to be a function of the previous sample mean $\nu_{t}(\MOtreatment)$, the sample mean of the current treatment $\bar{w}_{t+1}$, and an additional term to model the interaction of the dynamics and previous treatments $\nu_{t}(\MOtreatment) \bar{w}_{t}$. Consequently, this model is irrelevant of the variance $\rho_{t+1}(\MOtreatment)^2$. This is true when $g_t$ takes a simple nonlinear form $g_t(y_t, \vec{w}) = \alpha y_t + \beta w_{t+1} + \gamma y_t w_t$. We remark that FO-CMP essentially uses the first state evolution equation in (ref) and fails to extract informative signals from the second evolution equation.
HO-CMP (Higher-Order Causal Message-Passing): HO-CMP further introduces the second-order terms $(\bar{w}_{t+1})^2$ and $\hat{\rho}_t(\MOtreatment)^2$ to model the nonlinear treatment effects. It improves data efficiency by utilizing both state evolution equations. It also allows estimation of higher order terms in Taylor series of $g_t$.
While FO-CMP extends the estimation algorithm in shirani2024causal to accommodate experiments with more than two stages, HO-CMP introduces a new dimension to the estimation problem by incorporating second-order terms. This inclusion enhances data utilization, resulting in higher estimation efficiency in HO-CMP compared to FO-CMP.
\LinesNotNumbered \SetAlgoNoLine
In this section, we use synthetic experiments under simulated and real-world network interference patterns, to compare the performance of FO-CMP and HO-CMP estimators, outlined in Table (ref) and Algorithm (ref), with several benchmarks. First, we introduce the experimental design, benchmark estimators, interference patterns, and outcome specifications.
\paragraph{Experimental design.} We primarily focus on the staggered rollout design with $L$ distinct treated probabilities, denoted by $\pi^{(1)}, \cdots, \pi^{(L)}$, where $\pi^{(\ell)}$ increases monotonically with $\ell \in \{1,\ldots,L\}$. In the first $T^{(1)}$ periods, $\pi^{(1)}\times 100\%$ of units are in the treatment group. From $T^{(1)}$ to $T^{(2)}$ periods, $\pi^{(2)}\times 100\%$ of units are in the treatment group, and so forth. In the staggered rollout design, once a unit is allocated to treatment, it remains in the treatment group until the experiment concludes xiong2019optimal. In the appendix, we also consider the Bernoulli randomized design, where the treatment is re-randomized at every time period, allowing units to switch between the treatment and control groups throughout the experiment. We use two values of $T=40,\; 200$ and set $L=4$, with $(\pi^{(1)}, \pi^{(2)}, \pi^{(3)}, \pi^{(4)}) = (0.1, 0.2, 0.4, 0.5)$. In the appendix, we show the impact of increasing $L$ or the maximum treatment probability $\pi^{(L)}$.
\paragraph{Benchmark estimators.} We first present two benchmark estimators commonly used for treatment effect estimation, both in settings with and without network interference. The final estimator is designed specifically for settings with unknown network interference cortez2022staggered.
The first benchmark estimator is the standard difference-in-means (DM) estimator given by \[\widehat{\TTE}^\dm_{t} = \frac{\sum_{j=1}^N \Ooutcome{j}{t}{} \Otreatment{j}{t}{} }{\sum_{j=1}^N \Otreatment{j}{t}{} } - \frac{\sum_{j=1}^N \Ooutcome{j}{t}{} (1-\Otreatment{j}{t}{}) }{\sum_{j=1}^N (1-\Otreatment{j}{t}{}) } \, , \] which is the difference in average outcomes between treated and control units at each time period $t$.
The second benchmark is the standard horvitz1952generalization (HT) estimator given by \[\widehat{\TTE}^\hoth_{t} = \frac{1}{N} \sum_{j=1}^N \left[\frac{ \Ooutcome{j}{t}{} \Otreatment{j}{t}{} }{\pi_t} - \frac{ \Ooutcome{j}{t}{} (1-\Otreatment{j}{t}{}) }{1 - \pi_t }\right] \, , \] which weights observed outcomes by the inverse propensity score (i.e., $1/\pi_t$ or $1/(1-\pi_t)$).
The third benchmark estimator is the polynomial interpolation estimator (PolyFit) introduced by cortez2022staggered. PolyFit operates by obtaining estimates for the average of outcomes at equilibrium for $L$ treated probabilities $\pi^{(1)}, \cdots, \pi^{(L)}$, denoted by $\nueq(\pi^{(1)}),\ldots,\nueq(\pi^{(L)})$, then it utilizes Lagrange interpolation method and obtains a degree-$L$ polynomial approximation for the function $\nueq: [0,1]\to\mathbb{R}$ which can be used to estimate the equilibrium values under global control and treatment, $\hat{\nu}_{\mathrm{equil}}(0)$ and $\hat{\nu}_{\mathrm{equil}}(1)$. Finally, TTE is estimated by \[ \widehat{\TTE}^{\mathrm{polyfit}}_t = \hat{\nu}_{\mathrm{equil}}(1) - \hat{\nu}_{\mathrm{equil}}(0)\,. \] On the one hand, PolyFit does not need any knowledge of the interference network; however, it comes at the expense of having to grapple with two challenges. First, it may incur a high variance as $L$ increases due to fitting a high-degree polynomial. The second challenge is that it needs accurate estimates for each $\nueq(\pi^{(\ell)})$, which requires treating $\pi^{(\ell)}$ fraction of units for a long enough number of periods so that the outcomes reach an equilibrium. This can be achieved if the staggered roll-out design is performed over a long enough horizon $T$ with each $T^{(\ell)}$ sufficiently large, and then estimating each $\nueq(\pi^{(\ell)})$ by sample average of outcomes at time $T^{(\ell)}$. However, when such a lengthy experiment is not feasible, the estimates for $\nueq(\pi^{(\ell)})$ will be less accurate.
\paragraph{Interference networks.} We consider two networks (graphs). The first graph is a simulated random geometric graph model, studied by leung2022causal. The second graph is a social network of Twitch users rozemberczki2021twitch. In either scenario, we denote the adjacency matrix of the graph by $E \in \{0,1\}^{N \times N}$. For any $i$ and $j$, $E_{ij}$ equals $1$ if $j$ is a neighbor of $i$ and $0$ otherwise.
\paragraph{Outcome generating processes.} We consider two outcome specifications to model monotone and non-monotone interference patterns. Specifically, for both settings, we generate outcomes using the following specification: \[\outcome{i}{t+1}{} = \alpha + \beta \frac{\sum_{j=1}^N E_{ij} \outcome{j}{t}{} }{\sum_{j=1}^N E_{ij}} + \delta \cdot g\left( \frac{\sum_{j=1}^N E_{ij} \treatment{j}{t+1}{} }{\sum_{j=1}^N E_{ij}} \right) + \gamma \treatment{i}{t+1}{} + \epsilon^i_{t+1} \,, \] where in the first setting, $g(\cdot)$ is taken to be the identity function, i.e., $g(x) = x$ for any $x$. Therefore, $\outcome{i}{t+1}{}$ depends linearly on the fraction of treated neighbors, and we refer to this setting as the LinearInMeans outcome setting. This setting is widely studied in the causal inference literature cai2015social,eckles2016design,leung2022causal.
In the second setting, $g(\cdot)$ is specified by a periodic function, i.e., $g(x) = \sin(\uppi x)$ for any $x$. Therefore, $\outcome{i}{t+1}{}$, on average, first increases and then decreases with the fraction of treated neighbors, as visualized by the Ground Truth curve in panel (a) of Figure (ref). We refer to this setting as the Non-LinearInMeans outcome setting.
\paragraph{Results.} We compare FO-CMP and HO-CMP with the three benchmarks for estimating the TTE across the aforementioned outcome specifications and interference networks for long ($T=200$) and short ($T=40$) horizons. In each scenario, we perform 100 simulations of the synthetic experiment. The resulting distributions of ground truth and estimated TTEs are shown in Figures (ref)-(ref). All experiments were conducted on a MacBook Air with an Apple M1 chip and 16 GB of memory, with each setting taking about 15 minutes for 100 iterations. The key takeaways are as follows.
First, the DM and HT estimators exhibit significant bias across all cases. This is intuitive, as they estimate the TTE without accounting for the network interference.
Second, in the LinearInMeans outcome setting, FO-CMP and HO-CMP achieve low estimation error and minimal bias. This holds for both long experiment durations ($T=200$), where outcomes reach equilibrium, and short experiment durations ($T=40$), where outcomes have not yet reached equilibrium, as shown in Figures (ref) and (ref), respectively.
Third, as expected, PolyFit's dependence on accurate estimates for each $\nueq(\pi^{(\ell)})$ requires a large $T$ to reduce estimation bias. This is evident when comparing Figures (ref) and (ref): with a smaller $T$, PolyFit shows bias. This is also demonstrated in panel (b) of Figure (ref), where the red points—which represent sample averages of outcomes at $T^{(1)},\ldots,T^{(4)}$—have not yet converged and are slightly lower than their ground truth (equilibrium) values. This causes PolyFit's estimation of $\hat{\nu}_{\mathrm{equil}}(1)$ to be inaccurate, leading to a large bias. In contrast, HO-CMP, as shown in panel (c) of Figure (ref), is immune to this problem as it is designed to work with off-equilibrium data. Even with a larger $T$, the degree-$L$ polynomial estimation costs PolyFit with higher variance than both FO-CMP and HO-CMP, as shown in Figure (ref). Overall, this underscores the more efficient data utilization of FO-CMP and HO-CMP through their ability to leverage off-equilibrium data.
Fourth, in the Non-LinearInMeans outcome setting, HO-CMP achieves substantially lower estimation error compared to FO-CMP, as shown in Figures (ref) and (ref). This makes intuitive sense, as the higher-order terms in HO-CMP better capture the nonlinearity of $\nu_{\mathrm{equil}}(\pi)$ in $\pi$, while leveraging the additional data on sample variance dynamics over time, thereby enhancing the estimation accuracy.
Finally, the proposed estimation method demonstrates robustness across different experimental setups, including both LinearInMeans and Non-LinearInMeans outcome specifications. Additionally, robustness to graph structure---random versus Twitch graph---is evident from comparing the left and right plots in Figures (ref)-(ref). In Figure (ref) of Appendix (ref), we also demonstrate the robustness of the proposed methods to various parameters: the number of treatment probabilities $L$, the maximum treatment probability $\pi^{(L)}$, and the choice of experimental design (staggered rollout versus Bernoulli randomization).
Estimating causal effects under pervasive interference presents significant challenges sussman2017elements. Building on the causal message-passing framework of shirani2024causal, we incorporate higher-order moments of observed outcomes and treatment probabilities to estimate the total treatment effect, without requiring knowledge of the interference network. Our approach leverages machine learning techniques to extract informative patterns from these higher moments, enabling our estimator to capture complex counterfactual behaviors, including non-monotonic trends in outcome means relative to treatment proportions. While we demonstrate strong performance across various outcome specifications and network structures, the framework's applicability may be limited when multiple outcome observations are unavailable.