Higher-Order Causal Message Passing for Experimentation with Complex Interference

Authors: mohsen  bayati
yuwei  luo
william  overman
sadegh  shirani
ruoxuan  xiong
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 19:49:40

Abstract: Accurate estimation of treatment effects is essential for decision-making across various scientific fields. This task, however, becomes challenging in areas like social sciences and online marketplaces, where treating one experimental unit can influence outcomes for others through direct or indirect interactions. Such interference can lead to biased treatment effect estimates, particularly when the structure of these interactions is unknown. We address this challenge by introducing a new class of estimators based on causal message-passing, specifically designed for settings with pervasive, unknown interference. Our estimator draws on information from the sample mean and variance of unit outcomes and treatments over time, enabling efficient use of observed data to estimate the evolution of the system state. Concretely, we construct non-linear features from the moments of unit outcomes and treatments and then learn a function that maps these features to future mean and variance of unit outcomes. This allows for the estimation of the treatment effect over time. Extensive simulations across multiple domains, using synthetic and real network data, demonstrate the efficacy of our approach in estimating total treatment effect dynamics, even in cases where interference exhibits non-monotonic behavior in the probability of treatment.

Main paper: 
 introduction: 
 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) [Cox, 1958 , Rubin, 1978 , Manski, 1990 , Imbens and Rubin, 2015 , Sussman and Airoldi, 2017] .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 [Bond et al., 2012 , Blake and Coey, 2014 , Holtz et al., 2020 , Johari et al., 2022 , Bright et al., 2022] . 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 [Choi, 2017 , Cortez et al., 2022 , Li and Wager, 2022a] , 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 [Chen et al., 2024 , Agarwal et al., 2022 , Jia et al., 2024] , 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 [Manski, 2013 , Aronow and Samii, 2017 , Basse and Airoldi, 2018 , Karwa and Airoldi, 2018 , Forastiere et al., 2022] .Recently, Shirani and Bayati [2024] 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 [Mezard et al., 1986, Mezard and Montanari, 2009] and approximate message passing (AMP) [Donoho et al., 2009, Bayati and Montanari, 2011] 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 Shirani and Bayati [2024] , 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 Shirani and Bayati [2024] 's approach, which only employs the first moments for estimation. Second, while Shirani and Bayati [2024] 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 [Kohavi et al., 2020] .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 Cortez et al. [2022] , and a first-order CMP estimation, like the one by Shirani and Bayati [2024] . We emphasize that a large body of recent estimators, e.g., Jia et al. [2024] , 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.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 [Sobel, 2006 , Rosenbaum, 2007 , Hudgens and Halloran, 2012 , Tchetgen and VanderWeele, 2012 , Liu and Hudgens, 2014 , Kang and Imbens, 2016 , Viviano, 2020b , Bhattacharya et al., 2020 , Qu et al., 2021 , Auerbach and Tabord-Meehan, 2021 , Candogan et al., 2023 , Ugander and Yin, 2023] . When interference extends across clusters, standard estimators become biased. To address this, Eckles et al. [2016] 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 [Sussman and Airoldi, 2017] . Some recent studies combine the NIA with the availability of either a fully or partially observed interference structure [Leung, 2020 , Viviano, 2020a , Agarwal et al., 2022 , Belloni et al., 2022 , Li and Wager, 2022b] . Without prior knowledge of the interference structure, Cortez et al. [2022] consider low-degree polynomial interactions among units in the network. Leung [2022] 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 [Chin, 2018 , Jagadeesan et al., 2020 , Wang et al., 2020 , Li and Wager, 2022a , Agarwal et al., 2022 , Jagadeesan et al., 2020 , Leung, 2022] . 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 [Holtz et al., 2020 , Wager and Xu, 2021 , Munro et al., 2021 , Johari et al., 2022 , Harshaw et al., 2022 , Farias et al., 2022 , Bright et al., 2022 , Farias et al., 2023] . For example, Farias et al. [2022] 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 [Hudgens and Halloran, 2012 , Aronow and Samii, 2017 , Basse et al., 2019 , Jackson et al., 2020 , Sävje et al., 2021] . 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 [Li and Wager, 2022a, Boyarsky et al., 2023] . However, the problem of considering the dynamics of units' outcomes remains understudied [Arkhangelsky and Imbens, 2023] . 
setup and foundation: 
 Consider a system of N units indexed by i ∈ [N ] := {1, • • • , N } subject to a randomized experiment. The units are observed over a time horizon of T + 1 periods and for each t ∈ {0, 1, . . . , T }, we let W i t denote the treatment status of unit i during time period t. For simplicity, we consider a Bernoulli randomized design such that W i t ∼ Bernoulli(π t ). That is, at time t unit i receives the treatment with a probability of π t , corresponding to W i t = 1. Otherwise, unit i belongs to the control group and W i t = 0. In this context, we collectively define π = (π 0 , π 1 , . . . , π T ) as the experimental design. Then, following the potential outcome framework [Imbens and Rubin, 2015] , let Y i t (W) represent the potential outcome of unit i at time t, where W denotes the entire treatment allocation matrix, with W i t as the entry in row t and column i. Administering the treatment of unit i at time t according to w i t (as one realization of the random variable W i t ), we use w (as one realization of W) to show the matrix that captures the treatments of all units throughout the experiment; accordingly, we let y i t = Y i t (W = w) be the observed outcome of unit i at time t under the treatment assignment w:w =      w 1 0 w 2 0 . . . w N 0 w 1 1 w 2 1 . . . w N 1 . . . . . . . . . . . . w 1 T w 2 T . . . w N T      , y =      y 1 0 y 2 0 . . . y N 0 y 1 1 y 2 1 . . . y N 1 . . . . . . . . . . . . y 1 T y 2 T . . . y N T     .Observing (w, y), we are interested in estimating the TTE of the intervention, defined as below:EQUATIONwhere 1 and 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 [Jia et al., 2024 , Chen et al., 2024 , Viviano et al., 2023 , Yu et al., 2022 , Cortez et al., 2022] .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 [Holland, 1986] . Indeed, in Eq. ( 1), we can observe at most one of Y i t (1) or Y i t (0), and often, neither. 1 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, y and w, yielding accurate causal estimation under unknown network interference.In this section we provide a summary of the outcome specification and results of Shirani and Bayati [2024] that we utilize in the remaining. For t = 0, 1, . . . , T -1, we let g t : R × R T +1 → R be an unknown measurable function. We also use ⃗ W i = W i 0 , . . . , W i T ⊤ to denote the treatment assignment of unit i during the experiment. Accordingly, the treatment allocation matrix W is a T + 1 by N matrix with columns equal to ⃗ W i . Given potential outcomes Y j t (W) at time t and j ∈ [N ], their outcomes in time period t + 1 are specified byEQUATIONwhere G ij quantifies the impact of unit j on unit i at time t and ϵ i t is a zero-mean Gaussian noise with a variance of σ 2 e , accounting for measurement errors. In addition, we letG = [G ij ] i,j∈[N ]and refer to it as the interference matrix. Then, according to Eq. ( 2), the function g 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ν t (W) := lim N →∞ 1 N N i=1 Y i t (W), ρ t (W) 2 := lim N →∞ 1 N N i=1 Y i t (W) 2 -ν t (W) 2 . (3)Then, as shown by Shirani and Bayati [2024] , whenever the elements of the interference matrix G ij are i.i.d. Gaussian random variables with mean µ/N and variance σ 2 /N , under mild moment conditions on initial values Y i 0 , we have ν t+1 (W)EQUATION1 We can only observe y for one of exponentially many realizations of w.whereZ t ∼ N (0, 1) is independent from ⃗ W ∼ Bernoulli(π) (that is, W t ∼ Bernoulli(π t ) and ⃗ W = (W 0 , W 1 , . . . , W T ) ⊤) and the equalities hold almost surely. We note that the theory behind this result is rooted in the AMP literature, going back to Bolthausen [2014] , Bayati and Montanari [2011] . However, as Shirani and Bayati [2024] note, there is a major distinction between the AMP literature and the above setting: in the AMP literature, the matrix G is observed, and the aim is to construct proper functions g t for a completely different objective, which is studying the high-dimensional asymptotics of first-order algorithms. However, in the current context, the matrix G and functions g t are unknown and the goal is to estimate them.Considering Eq. ( 3), the equations in (4) 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 [Shirani and Bayati, 2024 ]. In the next section, we present an efficient algorithm to learn the state evolution dynamics outlined in Eq. ( 4). This method enables us to accurately estimate the TTE defined in Eq. ( 1) and its corresponding confidence interval.Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: We discuss the limitations of the work in Section 5. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
algorithm: 
 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 (4), thereby estimating counterfactual quantities while accounting for the impact of unknown network interference. To this end, by Eqs. ( 1) and (3), we rewrite the TTE as the difference of the sample means in the large limits:TTE t = ν t (1) -ν t (0).That means the problem of estimating the TTE is equivalent to estimating ν t (1) and ν t (0) using the observed data, denoted by (w, y). On the other hand, considering the state evolution equations in (4), the system state at time t + 1, denoted by (ν t+1 (w), ρ t+1 (w) 2 ), is a (nonlinear) function of the system state distribution at time t, characterized by (ν t (w), ρ t (w) 2 ) and ⃗ W , 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 (4) 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 (w, y). For this purpose, we fix the treatment assignment matrix w and defineνt (w) := 1 N N i=1 y i t , ρt (w) 2 := 1 N N i=1 y i t -νt (w) 2 , wt := 1 N N i=1 w i t , ⃗ w := ( w0 , . . . , wT ) ⊤ .In addition, let ⃗ ϕ = (ϕ k ) k∈[K] be a prespecified vector of measurable feature functions of current estimates of the sample mean νt (w), sample variance ρt (w) 2 , and the design w. We define x t to represent the feature vector as follows:x t = ⃗ ϕ νt (w), ρt (w), w := ϕ 1 νt (w), ρt (w), w , . . . , ϕ K νt (w), ρt (w), w ⊤ .Then, we formally propose learning the mapping f θ (•) defined by,EQUATIONWe summarize the method in Algorithm 1. Note in our experiment design we begin with all units under control by setting π 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 {π 1 , . . . , π T }. Specifically, in this paper, we consider the following estimators, as summarized in Table 1 .FO-CMP (First-Order Causal Message-Passing): This corresponds to the simple setting where ν t+1 (w) is assumed to be a function of the previous sample mean ν t (w), the sample mean of the current treatment wt+1 , and an additional term to model the interaction of the dynamics and previous treatments ν t (w) wt . Consequently, this model is irrelevant of the variance ρ t+1 (w) 2 . This is true when g t takes a simple nonlinear form g t (y t , ⃗ w) = αy t + βw t+1 + γy t w t . We remark that FO-CMP essentially uses the first state evolution equation in (4) 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 ( wt+1 ) 2 and ρt (w) 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 Shirani and Bayati [2024] 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.Algorithm 1: Higer-Order Causal Message Passing (HO-CMP)Data: Observed data (w, y), feature functions ⃗ ϕ = (ϕ k ) k∈[K] , machine learning model f θ (•) Step 1: Data processing for t ← 0 to T do νt (w) ← 1 N N i=1 y i t , ρt (w) 2 ← 1 N N i=1 (y i t -νt (w)) 2 , x t ← ⃗ ϕ νt (w), ρt (w) 2 , w end Step 2: Model Estimation Estimate f θ from data x t , νt+1 (w), ρt+1 (w) 2 t∈[T -1] , guided by (5). Step 3: Counterfactual Estimation ν0 (0) ← ν0 (w), ν0 (1) ← ν0 (w), ρ0 (0) 2 ← ρt (w) 2 , ρ0 (1) 2 ← ρt (w) 2 , TTE 0 ← 0 for t ← 0 to T -1 doCompute the features and predict the counterfactualsx t (0) ← ⃗ ϕ νt (0), ρt (0) 2 , 0 , x t (1) ← ⃗ ϕ νt (1), ρt (1) 2 , 1 νt+1 (0), ρt+1 (0) 2 ← f θ (x t (0)), νt+1 (1), ρt+1 (1) 2 ← f θ (x t (1)) Estimate the TTE TTE t+1 ← νt+1 (1) -νt+1 (0) end Result: TTE t t∈[T ]Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes]• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments: 
 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 1 and Algorithm 1, with several benchmarks. First, we introduce the experimental design, benchmark estimators, interference patterns, and outcome specifications.Experimental design. We primarily focus on the staggered rollout design with L distinct treated probabilities, denoted by π (1) , • • • , π (L) , where π (ℓ) increases monotonically with ℓ ∈ {1, . . . , L}.In the first T (1) periods, π (1) × 100% of units are in the treatment group. From T (1) to T (2) periods, π (2) × 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 [Xiong et al., 2024] . 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 (π (1) , π (2) , π (3) , π (4) ) = (0.1, 0.2, 0.4, 0.5). In the appendix, we show the impact of increasing L or the maximum treatment probability π (L) .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 [Cortez et al., 2022] .The first benchmark estimator is the standard difference-in-means (DM) estimator given byTTE dm t = N j=1 y j t w j t N j=1 w j t - N j=1 y j t (1 -w j t ) N j=1 (1 -w 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 Horvitz and Thompson [1952] (HT) estimator given byTTE ht t = 1 N N j=1 y j t w j t π t - y j t (1 -w j t ) 1 -π t ,which weights observed outcomes by the inverse propensity score (i.e., 1/π t or 1/(1 -π t )).The third benchmark estimator is the polynomial interpolation estimator (PolyFit) introduced by Cortez et al. [2022] . PolyFit operates by obtaining estimates for the average of outcomes at equilibrium for L treated probabilities π (1) , • • • , π (L) , denoted by ν equil (π (1) ), . . . , ν equil (π (L) ), then it utilizes Lagrange interpolation method and obtains a degree-L polynomial approximation for the function ν equil : [0, 1] → R which can be used to estimate the equilibrium values under global control and treatment, νequil (0) and νequil (1). Finally, TTE is estimated byTTE polyfit t = νequil (1) -ν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 ν equil (π (ℓ) ), which requires treating π (ℓ) 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 (ℓ) sufficiently large, and then estimating each ν equil (π (ℓ) ) by sample average of outcomes at time T (ℓ) . However, when such a lengthy experiment is not feasible, the estimates for ν equil (π (ℓ) ) will be less accurate.Interference networks. We consider two networks (graphs). The first graph is a simulated random geometric graph model, studied by Leung [2022] . The second graph is a social network of Twitch users [Rozemberczki and Sarkar, 2021] . In either scenario, we denote the adjacency matrix of the graph by E ∈ {0, 1} N ×N . For any i and j, E ij equals 1 if j is a neighbor of i and 0 otherwise. 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:Y i t+1 = α + β N j=1 E ij Y j t N j=1 E ij + δ • g N j=1 E ij W j t+1 N j=1 E ij + γW i t+1 + ϵ i t+1 ,where in the first setting, g(•) is taken to be the identity function, i.e., g(x) = x for any x. Therefore, Y 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 [Cai et al., 2015 , Eckles et al., 2016 , Leung, 2022] .In the second setting, g(•) is specified by a periodic function, i.e., g(x) = sin(πx) for any x. Therefore, Y 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 1 . We refer to this setting as the Non-LinearInMeans outcome setting.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 2 3 4 5 . 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 2 and 3 , respectively. Third, as expected, PolyFit's dependence on accurate estimates for each ν equil (π (ℓ) ) requires a large T to reduce estimation bias. This is evident when comparing Figures 2 and 3 : with a smaller T , PolyFit shows bias. This is also demonstrated in panel (b) of Figure 1 , where the red points-which represent sample averages of outcomes at T (1) , . . . , T (4) -have not yet converged and are slightly lower than their ground truth (equilibrium) values. This causes PolyFit's estimation of νequil (1) to be inaccurate, leading to a large bias. In contrast, HO-CMP, as shown in panel (c) of Figure 1 , 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 2 . 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 4 and 5 . This makes intuitive sense, as the higher- order terms in HO-CMP better capture the nonlinearity of ν equil (π) in π, 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 2 3 4 5 . In Figure 6 of Appendix A, we also demonstrate the robustness of the proposed methods to various parameters: the number of treatment probabilities L, the maximum treatment probability π (L) , and the choice of experimental design (staggered rollout versus Bernoulli randomization).  
conclusion: 
 Estimating causal effects under pervasive interference presents significant challenges [Sussman and Airoldi, 2017] . Building on the causal message-passing framework of Shirani and Bayati [2024] , 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. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: We provide the details of the computer resources in Section 4.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification:Guidelines: The research conducted in the paper adheres to the NeurIPS Code of Ethics.• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: We discuss broader impacts in Section 1.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: We cite the real network graph in Section 4.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 
 
 Appendix: 1. Claims Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? Answer: [Yes] Justification: Guidelines: Yes, the theoretical results in Sections 2 and 3 and the empirical results in Section 4 and Appendix A support the main claims made the abstract and introduction.• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.38th Conference on Neural Information Processing Systems (NeurIPS 2024). 


Decision: This work mirrors somewhat closely the work of Sadegh Shirani and Mohsen Bayati (https://arxiv.org/abs/2311.08340) and in light of that (beside concerns of plagiarism or self-plagiarism, since I don't see the authors' identities), the contribution of this submission remains a bit incremental.

Reviews: review: paper_summary: paper_summary: The authors motivate a regression framework using Causal Message Passing to estimate global average treatment effects under some unknown forms of interference. The simplest way to think about their proposed algorithm is that CMP is used to motivate particular sufficient statistics for the interference dynamics which are fed into a regression model to perform corrections and enable estimation.
strength_weakness: strength_weakness: The biggest problem with this paper is that it doesn't actually set out theoretical scope conditions for when it works. I'll be a bit more specific on two points:
1. The combination of the definition of TTE_t as the large sample limit combined with HO-CMP-I under Bernoulli randomization. Suppose W_t is determined by Bernoulli randomization (with a fixed probability, \pi across all time periods). Then as N \to \infty, then \bar{w}_t \to \pi. This means that the relationship between \bar{w}_t and \hat{\nu}_t will always be exactly zero. Even under the useful linear form discussed on line 228, it appears that this method could not possibly learn anything effectively. I suspect, in fact, that this would also be problematic in finite sample sizes (see below notes on evaluation). Critically, however, you haven't spelled out any of these scope conditions, so the paper does not make clear what is necessary to achieve success. And I will reiterate that this is using the very asymptotics gestured to in the paper and an experimental design which is explicitly discussed.
2. More specifically, I will return to line 228, which gestures at a generative model under which the proposed method (HO-CMP-I) should work well. But what does this mean? Should we expect the estimator to be unbiased? Should we expected it to be BLUE? Should we expect it to be semi-parametrically efficient? It would be ideal to have theoretical results around this, but barring that, these questions should be addressed directly and explicitly in the experiments.

On the subject of the experiments, I want to see more specifics than the charts you've provided give. Are the estimators unbiased? What are their MSEs? How do these evolve as a function of sample size? What are their convergence rates? In the absense of theoretical results, it is critical to have clear answers to these kinds of questions. 

Some minor points that I don't think matter much at all, but will share regardless:
- I think "higher-order" is an overly confusing way to define this approach. When I think higher-order, I think about something like a Taylor series: using more granular information about the setting to define a better estimator. In contrast, your approach is doing something more akin to higher-_aggregation_. I would consider using a term like "higher-level" or something which connotes this kind of aggregation.
- On line 211 you refer to "multi-label" ML models, but I think you mean "multi-task" models or "multivariate" models.
- What you've done here doesn't seem obviously like "machine learning." While I have no problem with linear regression, what it appears that you've done here
- line 264: you don't need to tell the reader what \pi is :)
questions: questions: see above
limitations: limitations: see above
paper_summary: paper_summary: This paper studies experimental interference, in which the treatment assignments of one unit affects the outcome of another. The majority of work in this area assumes interference acts through a network, and requires knowledge of that network to reduce bias in the resulting estimator. This paper takes a general approach in which interference occurs due to an unknown graph, and a unit's outcome can be affected by its neighbors' treatment status as well as their past outcomes. With such a flexible potential outcomes model, it is very difficult to estimate the total treatment effect in an unbiased way. The paper suggests an extension of the causal message passing algorithm that estimates the state evolution equations of the potential outcomes using higher-order functions of the past treatment status and outcomes, as opposed to the original CMP paper which used linear models of these quantities.
strength_weakness: strength_weakness: If I understand correctly, this work extends the CMP algorithm of Shirani and Bayati by using nonlinear estimators of the state dynamics. It is unclear to me that this delta is a significant enough improvement over the original CMP method to merit publication at Neurips.

I reviewed the Shirani and Bayati paper, and am concerned at how closely the writing, structure, and in some cases individual sentences of this paper parallel that work. I will leave it to the AC to make any recommendations about copyright / academic integrity, but here are the passages of most concern:
1. The phrase "Inspired by the literature on AMP algorithms, we refer to (4) as state evolution equations of the experiment" appears verbatim in the original CMP paper.
2. The related works of this paper appears to be a reworded and edited version of a subset of the related works of the CMP paper. In particular, the paragraphs beinning on lines 124 and 134 almost exactly mirror their counterparts in the CMP paper.
questions: questions: 1. Is the HO-CMP-I algorithm the same as that of Shirani and Bayati? If not, I would suggest including the original CMP algorithm as a baseline.
2. [Authors do not need to provide new experiments; discussion would suffice.] The HO-CMP algorithm outperforms polyfit as shown in Figure 1. It's unclear to me if this is due to increased regularization (maybe a 2-degree polynomial provides a better approximation to a half-sine than a 4-degree polynomial on the data in Figure 1) or an inherent property of the HO-CMP algorithm (maybe the claim is that not all types of interference are captured by the polyfit algorithm?). I would be curious if the authors had some intuition here, perhaps an example where polyfit wasn't obviously overfitting, but HO-CMP still outperformed polyfit.
limitations: limitations: Yes
paper_summary: paper_summary: This paper discusses the causal effect estimation under spatiotemporal interference characterized by a dynamic system. The authors are devoted to improving the causal message-passing framework by adding interaction terms of inputs in the summary function $f_\theta$ (a linear regressor). Synthetic and semi-synthetic experiments are conducted to show the performance of the proposed method.
strength_weakness: strength_weakness: 1. The comparison with the most related work [1] is very important but scarce, and I think the marginal contribution is little.  
The basic framework is all the same with [1], and the scanty innovation of this paper is the so-called "higher order", which is formulated as two more interaction terms and a square term. I don't think these additional input terms are necessary when given a non-linear model $f_\theta$, especially when the authors don't develop any theory that necessitates the simplicity of linear regression.  

2. The literature review is disordered, and there are apparent misunderstandings for certain important works. I list some examples in the following. 
- The taxonomy "Restrictions on interference structure" and "Treatment effect dynamics and temporal interference" overlap severely, wherein [2] studies the general temporal interference characterized by a MDP, while the authors list this paper with other two-sided marketplace paper in "Restrictions on interference structure".  
- I don't understand why [3] and [4] are discussed in the category of "Partial interference", since these two papers are important and don't belong to this category. 

3. More discussion on the basic setting is needed. I think neither the authors nor [1] discuss the concrete and applicable scenario of the dynamic system to show it's indeed meaningful, especially when they are agnostic to the network. I take network interference as an instance for the cross-sectional part. As the number of units increases, how the network topology change indeed impact the scenario, e.g. whether the degree of units increases with $n$ (Erdos-Renyi) or just keeps constant? The details behind the asymptotics indeed matter.



Ref. 
[1] Sadegh Shirani and Mohsen Bayati. Causal message passing: A method for experiments with unknown and general network interference. arXiv preprint arXiv:2311.08340, 2023.  
[2] Vivek Farias, Andrew Li, Tianyi Peng, and Andrew Zheng. Markovian interference in experiments. Advances in Neural Information Processing Systems, 35:535–549, 2022.  
[3] Christina Lee Yu, Edoardo M Airoldi, Christian Borgs, and Jennifer T Chayes. Estimating the total treatment effect in randomized experiments with unknown network structure. Proceedings of the National Academy of Sciences, 119(44):e2208975119, 2022.  
[4] Johan Ugander, Brian Karrer, Lars Backstrom, and Jon Kleinberg. Graph cluster randomization: Network exposure to multiple universes. In Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 329–337, 2013.
questions: questions: I have no questions.
limitations: limitations: No, the limitations of this paper are discussed inadequately, as discussed in weaknesses.
paper_summary: paper_summary: The paper considers the problem of total treatment effect estimation under interference, when the interference structure is unknown.  It proposes an estimation approach based on fitting a data-driven model for low-dimensional dynamics of the mean and variance of observed experimental outcomes over time, motivated by state evolution dynamics.  These state evolution equations come from approximate message-passing (AMP) applied to the potential outcome dynamics chosen in Equation 2.  To produce total treatment effect estimates, initial means and variances are observed from the experiment and propagated forward via the learned state evolution model.
strength_weakness: strength_weakness: - The authors should clarify the scope of their algorithm.  Equation 2 defines specific dynamics and does not meet the generality emphasized in the introduction and abstract around equilibrium and interference.  These dynamics do not obviously cover what may be occurring in real-world systems (e.g. lacking memory and being additive).  The authors motivate their paper as relaxing assumptions, but introduce their own alternative assumptions in Section 2.1.  The necessity of those assumptions is also undiscussed, leaving the range of applicability unclear.  For example, the interference matrix is assumed iid Gaussian but the experiments contain binary edges.  
- While past research Causal Message Passing is cited prominently in the paper, the paper could use more discussion over what is novel and different here.  To my understanding, the novel contribution is primarily the HO-CMP algorithms introduced in Section 3.  While the experiments demonstrate improvements by using HO-CMP, the idea of data-driven modeling of the state evolution equations appears a relatively straightforward extension.  
- The HO-CMP algorithms in Table 1 only leverage simple feature specifications that are low-dimensional. 
- Algorithm 1's computational complexity and scalability is unaddressed.
- Overall the paper could use increased clarity.  Many complex ideas are mentioned and occasionally described in insufficient detail.  For example, the authors do not justifying the necessity of considering state evolution instead of modeling Equation 2 directly.
questions: questions: 1.  Why only consider the mean of W_t as a feature in Table 1? 
2.  Equation 3: sums should be over i not n?
3.  What are the initial conditions of the counterfactual estimation part of Algorithm 1?
4.  In Equation 2, are the diagonal entries (i=j) treated the same as the off-diagonal entries?  What is the motivation for this?
limitations: limitations: Limitations are adequately addressed.

Comments: title: Global comment
comment: I thank the authors for their clarifications. I will retain my score.
Thank you for your response.
I thank the authors for the clarifications in the rebuttal, but my opinion on the paper overall is unchanged.
We thank the reviewer for carefully reading our paper and providing insightful comments. Below we provide responses to the weaknesses and questions you raised. 

Response to Weaknesses:

Thanks for helping us clarify the scope of our algorithm. In the introduction, we highlighted our algorithm can generalize the existing estimation methods under network interference in two aspects. First, our algorithm efficiently uses the data before the system reaches equilibrium. Second, our algorithm does not require knowledge of the network. Indeed Equation (2) encompasses these two aspects. We use data before equilibrium so that there are T periods in total (not just the last period at equilibrium). Moreover, we do not require the knowledge of the interference matrix G in Equation (2).
In addition, the properties of lacking memory and being additive can be covered in function g_t(). Specifically, lacking memory happens when the coefficient of Y_t^j is zerol and the coefficient of W_0^j, …, W_t^j are zerol. Being additive occurs when g_t() is a linear function of Y_t^j, W_0^j, …, W_t^j. 

The Gaussian assumption of the interference matrix is imposed to theoretically prove the state evolution in Equation (4). However, a key point of this paper is to show that the treatment effect estimation based on state evolution is accurate for general graphs, including the cases where the interference matrix does not follow the Gaussian distribution. Therefore, in the experiments, we consider the graphs with binary edges and illustrate the benefit of using our proposed method. 

Thanks for helping us be more clear about our contribution as compared to [SB23]. It is true that the main contribution is in Section 3 that provides more general and more accurate estimation methods of treatment effects. Our algorithm is more general in the sense that we allow for any number of \pi’s, whereas [SB23] only accounts for two different \pi’s in the randomized experiment. Our algorithm is more accurate because of two conceptual differences with respect to [SB23]. First, we use the state evolution of both the first and second moments (\nu_t and  \rho_t^2), while [SB23] only uses the state evolution of the first moment. Second, we use a data-driven approach to learn g_t(), while [SB23] considers a linear specification of g_t().
You are right that the feature specification in Table 1 only include first and second moments. We have explored the inclusion of higher-order moments, but the variance and MSE are generally higher than the three versions of HO-CMP in Table 1. 
Our algorithm runs in time O(NT). Once the mean and variance are calculated for every time t (that requires linear runtime in both N and T) in the step 1, the remaining steps only use the mean and variance. Therefore the algorithm is scalable with N (as in practice N is generally much larger than T).

Thanks for helping us be more clear. The state evolution equations are derived based on Equation (2) that hold for a broad class of unknown functions g_t() and unknown interference matrix G. As state evolution equations only involve summary statistics, they provide potential for developing efficient treatment effect estimation methods (especially, more efficient than [SB23]). This is exactly the problem studied in this paper. 

Response to Questions:

W_t^j is binary. The higher-order moments are therefore the same as the mean. 

You are right. We will correct it in the revision.

We set the initial conditions (i.e., \nu(0), \nu(1), \rho^2(0), and \rho^2(1) at time 0) as the mean and variance of observed outcomes at time 0 (i.e., \nu(0) and \nu(1) are set as the same; similarly, \rho^2(0) and \rho^2(1) are set as the same). We then compute \nu(0), \nu(1), \rho^2(0), and \rho^2(1) for time t=1, …, T. The difference between \nu(1) and \nu(0) is the estimated TTE at time t. As one can expect, the estimation of TTE is more accurate as t increases. 

You are right that the diagonal and off-diagonal entries are treated the same. It is possible to treat them differently and then separately identify the direct and spillover effects. However, as our goal is to identify TTE, it is sufficient to consider the model in Equation (2) and then the estimated effect will be the sum of direct and spillover effects.
Thank you for reading our work. As outlined in our global response, this research extends the estimation method of [SB23] by adopting a more data-driven approach. We will carefully revise the manuscript to ensure that the unique contributions of our work are clearly discussed and distinguished.

Regarding the literature review, we first emphasize that all the mentioned works have made significant contributions to the network interference problem. Specifically, [FLPZ22] is an excellent work focusing on experiments in dynamical systems. However, as clarified by the authors, their approach considers systems where treating some units impacts others through a limiting constraint, which is a **special case** of interference structure between units. By "Restrictions on interference structure," we refer to works that address the interference problem within a particular setting. Meanwhile, "Treatment effect dynamics and temporal interference" refers to studies that examine variations in the treatment effect over time. Additionally, [YABC22] and [UKBK13] are both crucial papers that we cited in the section related to "Partial interference," even though they do not adhere to the partial interference assumption. Indeed, these works effectively highlight the challenges associated with the broad applicability of studies in the "Partial interference" category and propose innovative solutions to mitigate these issues. We will ensure that the presentation is clearer to avoid these types of confusion.

About the core of the theoretical foundations of Causal-MP which proves, under certain assumptions, as N grows, sufficient statistics (means and variances) of the outcomes evolve according to state evolution equations. We also provide empirical simulations under network structures and outcome specifications that are studied in the prior literature with diverse network topology (e.g. Linear-in-means on random graphs studied by Leung and we generalize it to non-linear specifications and on real network data). For additional details please see Section 5 of [SB23].
Thank you for your thoughtful comments. Indeed, the reviewer is right, and as explained in our global response, this work builds on the foundation laid by [SB23] by extending their estimation method. We would also like to clarify that there was no intention to mirror the writing of [SB23]. However, it was necessary to discuss the problem setup and relevant literature. In the revision, we will ensure that the differences in the writing are clear and that the unique contributions of this work are highlighted.

To compare HO-CMP and PolyFit, it is important to note that HO-CMP utilizes more data by incorporating the evolution of the experiment over time. In contrast, PolyFit relies solely on the outcome observation at the equilibrium point where the treatment effect stabilizes. As a result of this additional data, HO-CMP can effectively mitigate the issue of overfitting that often affects PolyFit. Additionally, HO-CMP-i can be considered as a generalization of [SB23] estimator, allowing for more than two stages in the experiment.
1. We thank Reviewer p2A8 for helping us make this point more clear. The framework of HO-CMP enables us to learn an arbitrary functional form of \nu_{t+1}(w). In HO-CMP-i, we specifically model this as a linear function of \bar{w}_{t+1}. It is important to note that \pi varies during the experiments. For example, in the setting of Figure 2, we set (\pi^{(1)}, \pi^{(2)}, \pi^{(3)}, \pi^{(4)})=(0.1,0.2,0.4,0.5) for four intervals of duration 10, each. This approach allows us to learn the coefficient of \bar{w}_{t+1} in the linear function.

2. As we explained in the global response, the primary focus of this paper is to develop estimators and validate their performance in practical applications rather than provide theoretical analysis. The reason is that the theoretical treatment of valid questions (such as the ones raised by the reviewer) are highly challenging in this setting. We need to adjust our expectations as in the setting studied here (involving pervasive interference and unknown network structures) it is very difficult to estimate TTE. For example, take the Non-LinearInMeans outcome specification that we introduce as a generalization of the LinearInMeans (studied in the literature, e.g., Leung 2022). To the best of our knowledge, the proposed HO-CMP framework is the first method to efficiently utilize data for estimating TTE in this setting, and PolyFit and the estimator of [SB23] which is similar to HO-CMP-i being also the best benchmarks one could find (no other estimator to our knowledge can be applied here).

3. We thank Reviewer p2A8 for the comments. The expected convergence rate is O(1/\sqrt{N}), as discussed in [SB23]. We agree that an empirical demonstration of the rate would be a valuable addition to the paper that we plan to include in the revision.
We agree with the reviewer that the term “higher-order” might be confusing. It is intended to highlight that our function approximation of the dynamics incorporates additional terms. We will address all minor comments to clarify any potential ambiguities.
We thank the review team for reading our work. Here, we outline our responses by clarifying the contribution of the current work and defer more details to the reviewers' responses. This work builds on the foundation established by [SB23], extending their method in two directions by introducing the family of Higher-Order Causal Message Passing (HO-CMP) algorithms.

First, HO-CMP incorporates higher-order moments of unit outcomes, unlike [SB23]'s approach, which only employs the first moments for estimation. Second, while [SB23] focus solely on two-stage experiments, our work leverages the additional data provided by having more than two experimental stages, which aligns with the common practice in the tech industry of rolling out treatments through a sequence of experiments [KTX20]. 

Overall, this paper emphasizes a data-driven approach to the causal message passing framework and the efficient use of experimental data. Additionally, we introduce the 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.

[SB23] S. Shirani and M. Bayati, Causal Message Passing: A Method for Experiments with Unknown and General Network Interference 

[KTX20] Kohavi, R., Tang, D., and Xu, Y. (2020). Trustworthy online controlled experiments: A practical guide to a/b testing. Cambridge University Press.