Higher-Order Causal Message Passing for Experimentation with Complex Interference

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.