{'fraction_non_alphanumeric': 0.0372895134130763, 'fraction_numerical': 0.010416908605272326, 'mean_word_length': 4.677149261603375, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Heterogeneous Treatment and Spillover Effects under Clustered Network Interference

The bulk of causal inference studies rules out the presence of interference between units. However, in many real-world settings units are interconnected by social, physical or virtual ties and the effect of a treatment can spill from one unit to other connected individuals in the network. In these settings, interference should be taken into account to avoid biased estimates of the treatment effect, but it can also be leveraged to save resources and provide the intervention to a lower percentage of the population where the treatment is more effective and where the effect can spill over to other susceptible individuals. In fact, different people might respond differently not only to the treatment received but also to the treatment received by their network contacts. Understanding the heterogeneity of treatment and spillover effects can help policy-makers in the scale-up phase of the intervention, it can guide the design of targeting strategies with the ultimate goal of making the interventions more cost-effective, and it might even allow generalizing the level of treatment spillover effects in other populations. In this paper, we develop a machine learning method that makes use of tree-based algorithms and an Horvitz-Thompson estimator to assess the heterogeneity of treatment and spillover effects with respect to individual, neighborhood and network characteristics in the context of clustered network interference. We illustrate how the proposed binary tree methodology performs in a Monte Carlo simulation study. Additionally, we provide an application on a randomized experiment aimed at assessing the heterogeneous effects of information sessions on the uptake of a new weather insurance policy in rural China.


Motivation
According to Cox (1958), there is interference between different units if the outcome of one unit is affected by the treatment assignment of other units. In the case of policy interventions or socioeconomic programs, interference may arise due to social, physical or virtual interactions. For instance, in the case of infectious diseases, unprotected individuals can still benefit from public health measures undertaken in the rest of the population, such as vaccinations or preventive behaviors, because these interventions reduce the reservoir of infection (Nichol et al., 1995;Bridges et al., 2000), the vector of transmission (Binka et al., 1998;Howard et al., 2000) and the number of susceptible individuals (Broderick et al., 2008;Kelso et al., 2009;Kissler et al., 2020). In labor market, job placement assistance can affect job seekers using this service, but it can also have an effect on other job seekers who are competing on the same labor market (McKenzie and Puerto, 2015). In education, learning programs may spill over to untreated peers through knowledge transmission paths (Forastiere et al., 2019;Schiltz et al., 2019). In marketing, the exposure to an advertisement might directly affect the consuming behavior of the exposed individuals, and indirectly affect other individuals that are influenced by the consuming choices of people in their social network (Parshakov et al., 2020). If the exposure to the advertisement takes place in social media, the spillover effect might be also explained by the propagation of the advertisement to non-exposed users that are virtually connected to the targeted ones (Chae et al., 2017).
In the presence of interference, the effect of the treatment status of other units on one's outcome is usually referred to as spillover effect. In economics and social sciences there has been increasing interest in estimating spillover effects on networks in many different contexts: Duflo and Saez (2003) study the role of information and social interactions in retirement plan decisions in the academic community; Cai et al. (2015) assess the spillover effect of training sessions on uptake of weather insurance in rural areas of China; Muralidharan and Sundararaman (2015) study the aggregate effects of school choices while Imai et al. (2020) evaluate the effects of the Indian National health insurance.
Spillovers are a crucial component in understanding the full impact of an intervention at the population-level. In fact, the scale-up phase of a policy requires knowledge about the mechanism of spillover, and how this would take place in the population where the intervention will be rolled out. Information about spillovers of public policies can also support decisions about how best to deliver interventions, and could be used to guide public funds allocation. Indeed, the presence of beneficial spillover effects allows treating a lower percentage of the population, because the untreated individuals would still benefit from the treatment provided to the targeted sample. The use of spillover effects to save resources could be further improved if we were able to target specific individuals who would increase the overall impact of the intervention.
A targeting (or seeding) strategy aims at delivering the intervention to certain individuals such that the impact on the overall population is maximized (e.g., Valente, 2012;Kim et al., 2015;VanderWeele and Christakis, 2019;Montes et al., 2020). Typically, seeding strategies are designed in settings where either an element of the intervention (e.g., information, flyers or coupons provided during the intervention) or the outcome (e.g., the adoption of a behavior or a product) diffuse through the network. In these settings, the goal is the identification of the set of nodes in the network that, if targeted, would maximize contagion or diffusion cascades. To do so, seeding strategies are designed using information on the network structure and the dynamics of contagion or diffusion. This 'influence maximization' problem is NP-hard and computer scientist have developed approximate algorithms that usually rely on simplified contagion processes (Kempe et al., 2003).
Indeed, it is typically assumed that susceptibility to direct treatment and to others, as well as the influential power, are homogeneous across individuals.
Here we take a different perspective. First, we investigate spillover effects of a unit's treatment on other units' outcome without specifying the mechanism through which this might take place.
Second, we focus on the assessment of the heterogeneity of susceptibility to direct and indirect treatment. Understanding these heterogeneities can guide the design of targeting strategies aimed at making the interventions more cost-effective When spillovers are not present, this can be achieved by targeting those with the highest treatment effect. In fact, in the field of personalized medicine it is well known that individuals with different characteristics might respond differently to the treatment (e.g., Murphy, 2003;Chakraborty and Murphy, 2014;Kosorok and Laber, 2019). In the presence of interference, we also have that different people might be more or less susceptible to the treatment received by other units. This means that not only the treatment effect but also spillover effects are heterogeneous. An assessment of the heterogeneity of treatment and spillover effect is crucial not only for designing a cost-effective roll-out of the intervention within the targeted population, but it can also allow a generalization of its impact to other populations.

Related Work
In the past decade causal inference literature has experienced a growing interest in the mechanism of interference, leading to (i) the assessment of bias for causal effects estimated under the nointerference assumption (Sobel, 2006;Forastiere et al., 2020), (ii) the design of experiments to either avoid or assess interference (Angelucci and Di Maro, 2015;Baird et al., 2018;Kang and Imbens, 2016), and (iii) the estimation of casual effects under interference. Estimators for treatment and spillover effects have first been developed under the assumption of partial interference, allowing interference within groups but not across different groups (Hudgens and Halloran, 2008;Tchetgen and VanderWeele, 2012;Liu and Hudgens, 2014;Liu et al., 2016;Forastiere et al., 2016;Basse and Feller, 2018;Forastiere et al., 2019). However, the assumption of group-level interference is often invalid, too broad or not applicable. Hence, several works focus on the estimation of causal effects in the context of units interconnected on networks, both in randomized experiments (Aronow and Samii, 2017;Athey et al., 2018;Leung, 2020) and in observational studies (Ogburn et al., 2017;Sofrygin and van der Laan, 2017;Forastiere et al., 2018Forastiere et al., , 2020. In the context of social networks, even in randomized experiments where the treatment is randomized at the unit-level, exposure to other units' treatment is not. Therefore the propensity of being exposed to different levels of treatments among network contacts will depend on the network structure. Aronow and Samii (2017) developed a Horvitz-Thompson estimator to adjust for the imbalance in this propensity across units under different individual and contacts' treatment status.
In parallel to this field of research on interference, in recent year, thanks to the availability of increased computing power and large data sets, researchers have started to think about advanced data-driven methods to assess the heterogeneity of treatment effects with respect to large numbers of features. In this regard, there has been a large number of contributions in causal inference on subgroup analysis to investigate heterogeneous effect (e.g., Assmann et al., 2000;Cook et al., 2004;Rothwell, 2005;Su et al., 2009;Varadhan and Seeger, 2013;Ratkovic et al., 2017). However, the standard methods for subgroup analysis have several drawbacks. In particular, (1) they strongly rely on the subjective decisions on the specific variables defining the heterogeneous sub-populations; and (2) they fail to discover heterogeneities other than the ones that are a priori defined by the researchers. In addition, a data-driven method avoids potential problems related to cherry-picking the subgroups with extremely high/low treatment effects (Assmann et al., 2000;Cook et al., 2004).
Hence, many data-driven algorithms for the estimation of heterogeneous causal effects have been proposed in recent years (Hill, 2011;Foster et al., 2011;Su et al., 2012;Athey and Imbens, 2016;Wager and Athey, 2018;Starling et al., 2019;Lechner, 2019;Hahn et al., 2020). 1 The underlying aim of these methods is to detect 'causal' rules defining subsets of the study population where the treatment effect for that subgroup deviates from the average treatment effect. This is done by selecting the most important features and their values that define a partition of the covariate space (the tree) where the treatment effect is 'significantly' heterogeneous.
Among these algorithms, many rely on extensions of the standard Classification and Regression Trees (CART) algorithm (Friedman et al., 1984) and Random Forest (RF) algorithm (Breiman, 2001), and are adapted to different settings Imbens, 2015, 2016;Zhang et al., 2017;Lee et al., 2018;Wang et al., 2018;Bargagli Stoffi and Gnecco, 2018;Guber, 2018;Johnson et al., 2019;. In particular, in their seminal contribution Athey and Imbens (2016) develop the honest causal tree (HCT) methodology. HCT is a causal decision tree algorithm that is aimed at discovering the heterogeneity in causal effects through a single binary tree. HCT is constructed using a criterion function aimed at maximising the heterogeneity in causal effects at each split, while penalizing splits with higher variance in the estimated conditional effects. 2 In addition to a modification of the criterion function, Athey and Imbens (2016) introduce the concept of honest splitting. The authors propose to split the overall learning sample into a discovery (or training) set and an estimation set. The former is used to discover the heterogeneity in treatment effects, while the latter is used to estimate these effects in the discovered sub-populations. The different role of the two sets avoids potential spurious heterogeneity discovery due to overfitting of the learning sample. Building on the HCT, Wager and Athey (2018) and  introduce the causal forest (CF) and generalized random forest (GRF). Both CF and GRF are ensemble methods where multiple trees are built and combined to improve inference on the heterogeneous causal effects.
HCT and similar tree-based methodologies have already been applied to various scenarios for the discovery of heterogeneous effects of air pollution (Lee et al., 2018), employment incentives , job training programs (Cockx et al., 2019), development finance projects (Zhao et al., 2017), cardiovascular surgeries (Wang et al., 2017), cancer treatments (Zhang et al., 2017), and health insurance (Johnson et al., 2019). The wide usage of tree-based algorithms is due, in particular, to their ability to deal with non-parametric settings in an efficient and interpretable way. Indeed, these algorithms do not assume any specific shape of the treatment effect function. HCT and similar methodologies built on CART have been widely employed because of various attractive features: i.e., (i) they can deal with a large number of variables that are potentially responsible for the heterogeneity; (ii) they are simple to understand and visualize, easy to interpret, computationally scalable; and (iii) they can deal with non-linear relationships in the covariates without the need of data pre-processing. Nevertheless, tree-based methods for the estimation of heterogeneous causal effect have been developed ruling out the presence of spillover effects by assuming no-interference between units. On the other hand, the growing literature on spillover effects has focused on the estimation of population average spillover effects, neglecting potential heterogeneous spillover effects.
There have been few articles dealing with different types of heterogeneity in spillover effects. Forastiere et al. (2016Forastiere et al. ( , 2019 estimated the heterogeneity of spillover effects with respect to principal strata defined by the compliance behaviors in response to a cluster randomized treatment. However, the latent nature of these strata makes it difficult to effectively use the detected heterogeneity to design targeting strategies or to generalize the conclusions to a different population. Observed heterogeneity is instead studied in Arduini et al. (2014) and Arduini et al. (2019) where the focus is on the heterogeneity of peer effects from other units' outcomes across two specific groups, and the estimation relies on linear-in-means models and two-stage least squares. To the best of our knowledge there are no studies dealing with the heterogeneity of spillover effects on networks.

Contributions
In this paper, we bridge the gap between the aforementioned two bodies of causal inference literature by proposing a new algorithm for the discovery and estimation of heterogeneous treatment and spillover effects with respect to a large number of characteristics, including individual and network features. Our method is designed for randomized experiments affected by the presence of clustered network interference, that is, units are organized in a clustered structure, with no interactions between clusters and a network of connections within clusters (e.g., friendship networks within schools). In addition, randomization is assumed at the individual-level, resulting in treated and untreated units in the same cluster. Under this setting, spillover effects are confined within clusters and are assumed to take place on network interactions.
Our proposed method, network causal tree (NCT), builds upon the causal trees proposed by Athey and Imbens (2016), by modifying the splitting criterion to target treatment and spillover effects under interference. Splits are made so as to maximize the heterogeneity of the targeted causal effect(s), treatment and/or spillover effects, across the population. This criterion relies on the unbiasedness of the estimator of the effect(s) within each subset of the population. We first contribute to the existing literature on interference by developing an unbiased estimator for conditional treatment and spillover effects. We extend the Horvitz-Thompson estimator in Aronow and Samii (2017) to conditional causal effects under clustered network interference and we prove its consistency under the clustered network setting. This estimator is then used in our network causal trees to choose the binary splits that maximize the heterogeneity and to finally estimate the heterogeneous causal effects within the selected sub-populations.
In order to use the selected partition of the covariate space and the estimated treatment and spillover effects to guide policies, the heterogeneous sub-populations should be identified based on the causal effects that will be part of the decision rule. For instance, a policy that assigns the treatment to those who benefit the most from it requires the discovery of the subsets of the populations with heterogeneous treatment effect. Alternatively, a policy that is designed to target those who will respond to both their own treatment and the neighbors' treatment must identify sub-populations with high treatment and spillover effects. Hence, the proposed NCT methodology is optimized to detect the heterogeneity in treatment and spillover effects (i) either simultaneously, or (ii) separately. By reworking the criterion function of the seminal causal trees algorithm to account for interference, we allow the algorithm to detect heterogeneity in treatment and/or spillover effects. The discovery of the causal rules (the variables and their values defining heterogeneous subpopulations) representing the heterogeneity of one causal effect (treatment or spillover) is achieved by using a splitting criterion that maximizes the heterogeneity of that specific causal effect across the partition. On the contrary, if our goal is to identify a partition of the covariate space that can explain the heterogeneity of multiple causal effects, we propose the use of a composite splitting function that is designed to simultaneously maximize the heterogeneity in all the effects. This flexibility allows scholars and policy-makers to customize their investigations depending on their targeting goals. For instance, if a policy-makers wants to target individuals with the highest treatment effect and the lowest spillover effect (with the motivation that those with higher spillover effects can benefit from the treatment received by others), the NCT algorithm should be implemented with a composite splitting function aimed at detecting subsets of the population where both treatment and spillover effects are heterogeneous and the decision criterion can be applied. Conversely, if a targeting strategy is designed to target just individuals who would benefit the most from receiving the treatment, regardless of other people's assignment, a tree would be built using a single splitting criterion targeted to maximize the treatment effect heterogeneity. Similarly, a single criterion targeted to a spillover effect would be used in the case of targeting rules only involving that spillover effect.
It is important to note that the use of our algorithm to design implementation strategies is possible thanks to its high level of interpretability. Our network causal trees provide interpretable inference on heterogeneous treatment and spillover effects by discovering a set of causal rules that can be represented through a binary tree. As argued by Lee et al. (2018) and Lee et al. (2020) it is important to provide interpretable information on simple causal rules that can be targeted to improve policy effectiveness and to ensure that stakeholders and policy-makers understand (and, in turn, trust) the functionality of these models. Valdes et al. (2016) claim that a learning algorithm is interpretable if one can explain its classification by a conjunction of conditional statements (i.e., if-then rules). In this regard, tree-based algorithms based on if-then rules, such as the proposed NCT, are optimal for interpretability.
To assess the performance of the proposed NCT algorithm, we run a series of Monte Carlo simulations. In particular, we investigate the performance of the proposed algorithm with respect to two dimensions: its ability (i) to correctly identify the actual heterogeneous sub-populations and, (ii) to precisely estimate the conditional treatment and spillover effects. While the latter performance assessment is quite standard in the literature, the former is critical for interpretable algorithms for heterogeneous causal effects . Finally, we apply the proposed NCT algorithm to a randomized experiment conducted in China to assess the impact of information sessions on the purchase of a new weather insurance policy (Cai et al., 2015). Besides estimating the population average treatment and spillover effects (as already investigated in Cai et al. (2015)), our aim is to detect the strata of the population where one or both effects are heterogeneous and estimate these effects within these strata.
The remainder of the paper is organized as follows. In Section 2 we introduce the notation, setting, and assumptions that we employ throughout the paper. In Section 3 we define the conditional causal effects in a general partition of the covariate space and develop a Horvitz-Thompson estimator. Section 4 presents the proposed network causal tree algorithm, which is based on effect-specif or composite splitting functions for causal effects under interference. We then conduct a simulation study to assess the performance of the algorithm and estimator under different scenarios in Section 5 and we illustrate the application of the network causal tree on a randomized experiment in Section 6. Section 7 concludes the paper with a discussion of the proposed algorithm and directions for further research.

Notation and Setting
Let us consider a sample V of N units organized in K separate clusters. Let k ∈ K = [1, . . . , K] be the cluster indicator and let i = 1, . . . , n k be the unit indicator in each cluster k. Let now consider a connection structure such that units belonging to the same cluster might share a link whereas units belonging to different clusters are not connected. This network structure is represented by the graph G = (V, E), where V defines the set of nodes and E defined the set of edges, that is, the collection of links between each connected pair of nodes. A clustered network G is in turn an ensemble of K disjoint sub-graphs: G k = (V k , E k ), k = 1, . . . , K. The adjacency matrix A corresponding to the graph G, is a block-diagonal matrix with K blocks, A k , k : 1, . . . , K, where each element a ij,k is equal to 1 if there is a link between unit i and unit j in cluster k, that is, if the edge ij,k ∈ E k . Elements in A off the K blocks are equal to zero, indicating no links between units belonging to different clusters.
Let now W ik ∈ {0, 1} be a binary variable representing the treatment assigned to unit i in cluster k and let Y ik be the observed outcome. We denote by W k and Y k the treatment and outcome vectors in each cluster k. Similarly, W and Y denote the treatment and outcome vectors in the whole sample. It is worth noting that the treatment status is not the same for all units belonging to the same cluster. This treatment allocation corresponds to a unit-level randomization where treatment is assigned to each unit of a cluster independently (as in a Bernoulli trial) or with some level of dependency (as in a completely randomized trial) based on a probability distribution P (W k ) (see Section 2.3 for further details). Moreover, for each unit ik we observe a vector X ik of P covariates (or pre-treatment variables) that are not influenced by the treatment assignment. The vector of covariates might include individual characteristics (e.g., age, sex, socio-economic status, ...), cluster-level characteristics (e.g., cluster size, location, ...), as well as network characteristics representing aggregated individual characteristics (e.g., average age or proportion of males and females, ...) or the network topology (e.g., degree, centrality, transitivity, ...).

Clustered Network Interference
Following the potential outcome framework (Rubin, 1974;Holland, 1986), we denote by Y ik (w) the potential outcome that unit i in cluster k would experience if the treatment vector W in the whole sample were w, with w ∈ {0, 1} N . Under the assumption of no-interference, the potential outcome could be indexed only by the individual treatment assignment W ik , that is, Y ik (W ik = w).
In combination with the assumption of consistency, this assumption is known as Stable Units Treatment Value Assumption (SUTVA) (Rubin, 1986). The no-interference assumption is clearly violated in many real-world scenarios. For instance, the evaluation of the effect of vaccines should take into account the fact that unprotected individuals would also benefit from high vaccination coverage around them.
Here, we focus on a particular type of interference: clustered network interference. As we will show later, focusing on this type of interference is critical to ensure asymptotic properties of the estimator for conditional causal effects as well as to allow the network causal tree to divide the sample into a training (or discovery) set and an estimation set (and a testing set, if applicable). 3 The assumption of clustered network interference implies that: i) interference is restricted to nodes of the same cluster and interference between clusters is ruled out, that is, one's outcome is only affected by the treatment received by units belonging to the same cluster; ii) one's outcome is affected by a weighted function of the treatment status of potentially all units in their own cluster, with weights depending on the presence and possibly the value of network links.
Let W k/i be the vector collecting the treatment status of all units in cluster k except unit i.
Let g(·) : {0, 1} n k −1 −→ ∆ ik be a function that maps a cluster assignment vector W k/i to an exposure value. Without loss of generality, we define it as a function of the dot product between the cluster assignment vector and a vector of weights δ i (A k , X k ), which in turn depends on the adjacency matrix A k and the covariate matrix X k , i.e., g(W k/i , δ i (A k , X k )) = f (W k/i ·δ i (A k , X k )).
For instance, the function g(·) could result in the number or proportion of treated units in a cluster.
In this case the weight vector would be equal to δ i (A k , X k ) = 1 n k −1 or δ i (A k , X k ) = ( 1 n k −1 ) n k −1 , respectively. Alternatively we could use the adjacency matrix to compute the geodesic distance d(i, j) between each pair of nodes in cluster k and let g( is similar to the 'effective treatments' function in Manski (2013) and the 'exposure mapping' function in Aronow and Samii (2017), although it applies to the cluster treatment vector only. To ease notation, throughout we will omit the weight vector δ i (A k , X k ) in the function g(·).
We can now formalize the clustered network interference assumption as follows.
Assumption 1 states that the outcome of a unit i in cluster k depends on the individual treatment W ik and a function of the treatment status of the other members of cluster k, i.e., g(W k/i ), regardless of the specific treatment status of each member. This assumption can be viewed as an intermediate assumption between (i) assuming no interference and (ii) making no assumptions about the nature of interference. In a way, it is similar to the partial interference or the stratified interference in Hudgens and Halloran (2008), which are special cases of the clustered network interference assumption, with g(W k/i ) = W k/i and g(W k/i ) = n k j=1 W ij . Let G ik = g(W k/i ), referred to throughout as network exposure. Under Assumption 1, each unit has |∆ ik | × 2 potential outcomes, which we can write in terms of the individual treatment and the network exposure exposure as Y ik (w, g), representing the potential outcome of unit ik under We also assume the following consistency assumption: This assumption rules out different versions of the treatment and different ways in which a value of the network exposure can affect the outcome of a particular unit. Under a 'finite sample perspective', we assume the potential outcomes of each unit to be fixed but unknown, except for the observed Y ik (W ik , G ik ). Therefore, the only source of randomness in the potential outcomes is given by the random assignment to the treatment and the random network exposure induced by the random cluster assignment.
Assumptions 1 and 2 together are alternative to SUTVA when interference is present and is limited to within clusters. When the weight function δ i (A k , X k ) is such that elements δ ij (A k , X k ) = 0 if j ∈ V k : a ij,k = 0, that is, units that are not directly connected to unit i receive a weight equal to zero, then interference is limited to the neighborhood N ik of each unit, with N ik = {j ∈ V k : a ij,k = 1}. In this case, Assumptions 1 and 2 correspond to the SUTNVA Assumption in Forastiere et al. (2020). We denote by N g ik the set of units defining the network exposure, that is, In most of the literature on spillover effects this set is either the cluster k (Hudgens and Halloran, 2008) or the neighborhood of unit i (Forastiere et al., 2020). Alternative specifications are also possible and might involve higher-order neighbors.
For the purpose of assessing effect heterogeneity using tree-based methods, we will further make the following assumption: Assumption 3 (Discrete Network Exposure). We consider a discrete exposure mapping function Z is the set of integers. This assumption implies that the network exposure G ik is a discrete variable For instance, we can define a binary network exposure based on a threshold function applied to the number of treated neighbors: where q is a threshold. Hence, the g(·) exposure mapping function behaves like a threshold function which sums the elements of its argument (i.e the treatment assignment vector that characterizes the neighborhood of each unit) taking the value of 1 if the resulting value exceeds a certain threshold (e.g., at least one treated neighbor is treated, the majority of the neighbors are treated, ...). In our simulation study as well as in the application we have chosen the following definition: G ik = 1( j∈N ik W jk ) ≥ 1 , that is, the network exposure is 1 if at least one network neighbor is treated. As a consequence, both the individual treatment and the network exposure are defined as binary variables, W ik ∈ {0, 1} and G ik ∈ {0, 1}. It follows that the support of the joint treatment variable (W ik , G ik ) is finite and comprises four possible realizations, given by the combination of the two A discrete network exposure is crucial for our causal tree algorithm, at least in the version proposed in this paper. Indeed, the algorithm relies on the presence of enough observations for each treatment and exposure value to allow the estimation of the causal effects. Depending on the stopping rule which might rely on the accuracy of the estimation of conditional effects or on the number of observations (see Section 4), if the sample size is not large enough with respect to the number of categories of the network exposure and/or its distribution is non-uniform and highly skewed, the network causal tree algorithm might result in a tree with low depth and low granularity, that is, with highly heterogeneous causal effects even within the terminal leaves. Therefore, the maximum number of categories for the network exposure depend on the sample size, the number of covariates and their nature, as well as on the extent of the heterogeneity in the causal effects.

Unit-Level Randomization and Induced Marginal and Joint Distributions
In this work, we consider an experimental design with a unit-level randomization of the treatment, which is independent between clusters but might be dependent within them. Therefore, the treatment vector W is a random vector with probability distribution P (W = w) and the following assumption holds.
where W k is the treatment vector in each cluster k.
We denote by π W ik the unit-level probability that W ik is equal to 1, under the experimental design in place. In a randomized experiment π W ik is known. In the case of a Bernoulli trial, where each unit is independently assigned to the individual treatment, π W ik is constant and equal to α. 4 An example of a design with randomization independent between clusters but dependent within clusters is that of a completely randomized experiment taking place in each cluster. In this case, π W ik would be equal to m k /n k , where m k is the fixed number of treated units, and the treatment assignment for each unit does depend on the treatment status of other units.
Since the network exposure is a deterministic function g(·) of the cluster assignment vector W k/i , then the randomization distribution P (W = w) induces, together with the definition of the function g(·), a probability distribution of the vector of network exposures G in the whole sample. Hence, the probability for a unit of being exposed to a specif value of the network exposure G ik = g given the individual treatment w, denoted by π G|W ik (g|w), is known and can in principle be computed from the probability distribution P (W). Note that, we can drop the dependency from the individual treatment and write π G ik (g) when the randomization is independent between units.
Let ∆ = {0, 1} × ik∈V ∆ ik be the domain of the joint individual and network treatment status, that is, (w, g) ∈ ∆. Let π ik (w, g) denote the marginal probability for unit ik of being assigned to individual treatment w and being exposed to the network status g. This is equal to the expected proportion of assignment vectors w inducing an individual treatment w and a network exposure g: This marginal probability is a crucial component of the Horvitz-Thompson estimator for causal effects under network interference. For instance, if the experimental design is a Bernoulli trial with unit-level probability α and the network exposure is defined by a threshold function on the neighborhood as in Equation 1, then the joint probability could be computed as follows: where N ik is the number of neighbors ('degree') of unit ik.
To deal with well-defined potential outcomes, we must assume that each unit has a non-zero probability of being exposed to each (w, g): Assumption 5 (Positivity). π ik (w, g) > 0 ∀i ∈ N , k ∈ K and ∀(w, g) ∈ ∆.
When π ik (w, g) = 0 for some units, then the average potential outcomes and causal effects involving these values z and g must be restricted to the subset of units for whom π ik (w, g) > 0. For instance, if the network exposure is defined as in Equation 1, then the positivity assumption is violated for units who cannot be exposed to a value g, that is, those with a degree N ik lower than the threshold q. Consequently, the analysis must be restricted only to the subset of the population satisfying the positivity criterion.
The estimator that we propose below also requires the so-called pairwise exposure probabilities, which describes the joint probability for pairs of units of being exposed to a given individual treatment and network status. Hence, given specific exposure conditions (w, g) and (w , g ) , a pairwise exposure probability, denote by π ikjh (w, g; w , g ), quantifies the probability that the two events (W ik = w, G ik = g) and (W jh = w , G jh = g ) occur, i.e., π ikjh (w, g; w , g ) = P (W ik = w, G ik = g, W jh = w , G jh = g ). In general, this can be written as: Under the event of both units being exposed to the same condition (w, g) we denote the pairwise exposure probability by π ikjh (w, g).
In the case of an experimental design assigning treatment independently between clusters, under the clustered network interference the two events (W ik = w, G ik = g) and (W jh = w , G jh = g ), with k = h, are independent and the pairwise exposure probability equals the product of the two joint probabilities: π ikjh (w, g; w , g ) = π ik (w, g) × π jh (w , g ). In the case of a Bernoulli trial and the network exposure defined on the neighborhood only, this is also true for units i and j belonging to the same cluster, i.e., k = h, but are not connected and do not share any neighbors. In general, Even under independent treatment assignment, if N wg ik ∩ N wg jk = 0 the joint treatment of the units ik and jk will be dependent and π ikjk (w, g; w , g ) = π ik (w, g)π jk (w , g ) = 0. Note that if π ik (w, g) or π jh (w , g ) = 0 ⇒ π ikjh (w, g; w , g ) = 0, but not the reverse. Indeed, the joint probability of the two events (W ik = w, G ik = g) and (W jh = w , G jh = g ) might be zero if the network exposures G ik and G jh are defined on two subsets of units that coincide or include jh and ik, respectively. For example, if the network exposure is defined as in Equation 1 with threshold equal to 1 (i.e., having at least one treated neighbor) then, if unit ik is treated, with W ik and belongs to the neighborhood N jh of unit jh, the network exposure G jh cannot be 0.

Conditional Treatment and Spillover Effects
Our ultimate goal is to detect the regions of the covariate space exhibiting a high level of heterogeneity in the causal effects and estimate the causal effects of interest in these heterogeneous regions. In this section, we will focus on the definition and estimation of conditional treatment and spillover effects and we will assume that the heterogeneous regions that we want to investigate have already been identified, either a priori according to subject-matter knowledge or thanks to data-driven methods.
Let us denote with Π a partition of the covariate space X into M non-overlapping regions: m=1 m = X , and with (x, Π) : X → Π a function that maps each vector x of the covariate space into a region. Let N ( m ) be the size of each region m , with m = 1, . . . , M , and let V k ( m ) be the subset of units belonging to region m in cluster k, with k = 1, . . . , K. In the machine learning literature on CART, these non-overlapping regions are referred to as leaves.
For consistency, throughout we will use this terminology, regardless of whether the partition Π has been a priori defined or is the result of a tree-based algorithm. In addition, to ease notation, we will drop the reference to the partition Π from the mapping function (·).
When units are organized in a network, it is worth noting that a partition Π of the covariate space divides the sample units into sub-populations according to similarities in their characteristics, regardless of their network distance. Hence, two connected units might belong to different regions of the partition. However, in an homophilous network, where the probability of forming a link depends on the similarity in certain features and, hence, connected units are likely to share similar characteristics, a partition of the the covariate space is also likely to cluster connected units together. Given a partition Π, we now define conditional average potential outcomes under each individual treatment and network exposure condition (w, g) ∈ ∆. For the subset of units S m with covariate vectors x ∈ X that are mapped to the same region by the function (x), i.e., S m = {ik ∈ V : (X ik ) = m }, we define the leaf-specific average potential outcome under treatment and exposure condition (w, g) ∈ ∆ as follows: Note that µ (w,g) ( (x)) is a sample average, that is, it is the average potential outcome for all units in the sample V with a covariate vector mapped to the same region (x).
Leaf-specific conditional average causal effect (CACE) can be defined by comparing average potential outcomes under two different conditions: We denote by T the set of possible contrast we are interested in. For instance if both the individual treatment and the network exposure are binary, then T ⊆ {(1, 0, 0, 0), (1, 1, 0, 1), (0, 1, 0, 0), (1, 1, 1, 0), (1, 1, 0, 0)}. We define as leaf-specific treatment effects causal contrasts τ (w,g;w ,g) ( (x)) that keep the network exposure fixed at a level g while changing the individual treatment from w to w, that is, when g = g . These represent causal effects of receiving the treatment while the treatment status of all other units is kept fixed or is mapped to the same network exposure g. On the contrary, we define as leaf-specific spillover effects causal contrasts τ (w,g;w,g ) ( (x)) that keep the individual treatment fixed at a level w while changing the network exposure from g to g, that is, when w = w . These spillover effects can be seen as causal effects of a change in the treatment status of other units such that the network exposure also changes, while the individual treatment status is kept fixed. It should be emphasized that µ (w,g) ( (x)) corresponds to a unit-level intervention setting the treatment and network exposure of each unit to specific values. The focus on these type of average potential outcomes, as opposed to the ones based on population-level hypothetical interventions as in Hudgens and Halloran (2008), is due to our purpose of investigating heterogeneous responses to the individual treatment and network status across units with different characteristics. If one is interested in assessing the heterogeneity of the average response to the network exposure resulting from a hypothetical treatment allocation, our approach could be extended to marginalized causal effects as the ones in Forastiere et al. (2020).

Estimator for leaf-specific CACE
Here we develop an Horvitz-Thompson estimator for leaf-specific conditional average causal effects.
The derivation of the proposed estimator builds upon the estimator for average causal effects under network interference proposed by Aronow and Samii (2017).
Following Horvitz and Thompson (1952) and Aronow and Samii (2017), a design-based estimator for the leaf-specific average potential outcome under individual treatment w and network exposure g, µ (w,g) ( (x)), can be expressed as: where π ik (w, g) denotes the probability of a given unit ik, that belongs to the leaf (x) (in the partition Π), to be exposed to the treatment condition (w, g).
The variance estimator of µ (w,g) ( (x) can be expressed as: This expression extends the variance estimator derived in Aronow and Samii (2017) (Equation 7) to the case of conditional average potential outcomes and clustered network interference. In fact, the second term in (8) includes the covariance between the individual treatment and network exposure of two units belonging to the same leaf (x). Under an experimental design with independent treatment allocation between clusters and under clustered interference such covariance between two units belonging to different clusters is zero and the second term should be restricted to units j in the same cluster as i. In addition, the covariance between the joint treatment of two units is non-zero if the set of units defining the network exposure -e.g., the whole cluster or the unit's neighborhood -is shared between them or includes them. Formally, if N wg ik ∩ N wg jk = 0 the joint treatment of the units ik and jk will be dependent, that is, π ikjk (w, g) − π ik (w, g)π jk (w, g) = 0 .
Hence, two units belonging to the same leaf are more likely to have intersecting sets N wg ik and N wg jk (e.g., shared neighbors) if the sets are homogeneous, that is, units belonging to these sets share similar characteristics. Settings with homophilous networks are investigated in the appendix.
An estimator for the leaf-specific conditional average causal effect of the exposure condition (w, g) compared to the configuration (w , g ) can be written as: The estimated variance of the estimator τ (w,g;w ,g ) ( (x)) can be decomposed as follows: with the covariance estimator taking the following expression for the case when π ikjk (w, g; w , g ) > 0 ∀i, j, k: Further details about the variance estimator of leaf-specific CACE can be found in Appendix B.

Properties of the Horvitz-Thompson Estimator
Here we will describe the properties of the Horvitz-Thompson estimator of leaf-specific causal effects. Asymptotic results will rely on a growth process that is commonly assumed with cluster data. In particular, we consider a sequence of nested samples V of size N, where V consists of K separate clusters V k of size n k , k = 1, . . . , K. We let the sample size N −→ ∞ by letting the number of clusters go to infinity, i.e., K −→ ∞, while the cluster size n k , k = 1, . . . , K remains fixed.
Proof. Proof in Appendix A.
The unbiasedness of the estimator of leaf-specific CACE is conditional on the partition Π and the function (·). When building causal trees to assess the heterogeneity of causal effects, we will rely on this property to derive the splitting criterion and to estimate leaf-specific causal effects. 5 The proof follows directly from the unbiasedness of the Horvitz-Thompson estimator. A conservative estimator for the case when π ikjk (w, g) = 0 for some units can be found in Appendix B.
A proof follows from Aronow and Samii (2017). The restriction on the covariances does not change the proof.

Proof. See Appendix A.
Note that cluster network interference and independent treatment allocation between clusters ensure that the amount of dependence across units is limited. This limited independence is the condition required to ensure consistency (Aronow and Samii, 2017). 6 Proposition 5 (Asymptotic Normality). Given an cluster independent design and the clustered network interference assumption, then: ) .
An independent treatment allocation between clusters and the clustered network interference ensure the limited dependence condition required in Aronow and Samii (2017). This condition allows us to rely on a central limit theorem derived via Stein's method (Chen and Shao, 2004) to achieve the asymptotic normality of the estimator. The variance estimators will depend on the size of the sample belonging to leaf (x) in each cluster, i.e., n k ( (x)) ≤ B( (x)) ≤ B, k = 1 = . . . , K, and the maximum conditional degree D( ) = max ik∈V : ). Given that these quantities are bounded, we can show that V V τ (w,g;w ,g ) ( (x)) = O(1/N ( (x))), that is, the rate of convergence will be 1/ N ( (x)), with N ( (x)) ≤ KB( (x)) (the proof follows the one in Aronow and Samii (2017)).
It is worth noting that, thanks to the asymptotic growth assumed here, Proposition 5 would still hold if the covariate vector X included network covariates defined both both at the cluster-level or at the individual-level. This result would allow us to investigate the heterogeneity of causal effects with respect to network characteristics, including variables defining a cluster network structure or the structure of the network neighborhood around a node.

Network Causal Trees for Heterogeneous Causal Effects under Clustered Network Interference
In the previous section we have introduced and developed an estimator for causal effects conditional on sub-populations of units defined by a partition Π of the covariate space X . Here we develop a data-driven machine learning algorithm to identify the partition Π aimed at investigating the heterogeneity in the effects of interest.
Our proposed algorithm, named Network Causal Tree (NCT), builds upon the Causal Tree (CT) algorithm introduced by Athey and Imbens (2016), which in turn finds its roots in the Classification and Regression Tree (CART) algorithm (Friedman et al., 1984). CART is a widely used nonparametric method to partition the feature space. It relies on a tree-based algorithm which recursively splits the sample. In particular, trees are constructed by recursively partitioning the observations from the root (that contains all the observations in the learning samples) into two child nodes. This procedure is repeated until the tree reaches the final nodes which are called leaves. Because each node is always split into two sub-nodes, these trees are called binary trees.
Binary trees are called regression trees when the outcome is a continuous variable, while they are called classification trees when the outcome is either a discrete or a binary variable. The aim of the tree construction is to identify heterogeneities in the relationship between the observed outcome and the features to best predict the outcome variable. Therefore, splits are made with the aim of minimizing the prediction error. With this aim different splitting criteria could be specified. For additional details on CART, we refer to the seminal paper by Friedman et al. (1984). Figure 3 illustrates an example of binary partitioning in a simple case with just two predictors x 1 ∈ [0, 1] and x 2 ∈ [0, 1].
x 1 < 0.6 x 2 > 0. Building on CART, Athey and Imbens (2016) developed a causal decision tree algorithm with the aim of detecting the heterogeneity of causal effects. In particular, they modified the splitting function to minimize the estimation error of conditional effects. Moreover, Athey and Imbens (2016) introduced honest inference by using a sub-sample to build the tree (training or discovery sample) and a separate sub-sample to perform inference (estimation sample). This sample-splitting approach is transparent and efficient even in high-dimensional settings.
Our proposed NCT differs form the standard causal tree algorithm in two critical aspects: (i) it estimates heterogeneous causal effects -both treatment and spillover effects -in the presence of clustered network interference, and (ii) it possibly models heterogeneity with respect to more than one effect at the same time through a composite splitting criterion. In this Section, we describe and motivate the splitting criteria for our NCT algorithm (Subsection 4.1) and its detailed structure (Subsection 4.2).

Splitting Criteria
The NCT algorithm is built to detect and estimate heterogeneous treatment and spillover effects, in the presence of clustered network interference. Moreover, NCT is able to discover the heterogeneity with respect to more than one estimand. Here, we present the three criteria that rule the splitting procedure of NCT, targeted to single effects or multiple effects.

Single splitting criteria
Let P be the space of partitions. Given a causal effect τ (w,g;w ,g ) , we can use recursive splitting to look for the best partition Π ∈ P with respect to a splitting criterion Q (w,g;w ,g ) (Π). Formally: Given our goal of describing the relationship between the causal effect and the covariate space and detecting subsets that exhibit a high level of heterogeneity, we can define a splitting criterion that maximizes accuracy in the prediction of conditional effects τ (w,g;w ,g ) (X ik ) in the whole sample V.
This translates into the minimization of expected value of the mean square error (MSE): where the expected value is taken over the sampling distribution. When this splitting criterion is used to select the partition Π, we maximize the function in (13) evaluated in the sample used to build the tree, i.e., the training set. In this case, in the machine learning literature the objective function is referred to as the in-sample splitting function and we denote this by Q in (w,g;w ,g ) (Π). As opposed to the EMSE of the observed outcome prediction, the true causal effect τ w,g;w ,g (x) is unknown. However, we can use the training data to estimate the EMSE for the in-sample spitting rule. Thanks to the unbiasedness of the estimator τ (w,g;w ,g ) ( (x), Π) with respect to the population causal effect E[τ (w,g;w ,g ) (X ik )|X ik ∈ (x)] (Proposition 6 in Appendix A), following Athey and Imbens (2016) we can estimate the EMSE as follows: where K tr is the subset of clusters belonging to the training set and N tr is the sample size. 7 Therefore, the maximization of this splitting function results in the maximization of the heterogeneity across leaves. In fact, if two sub-populations 1 and 2 have a different causal effect τ (w,g;w ,g ) , i.e., τ (w,g;w ,g ) ( 1 ) = τ (w,g;w ,g ) ( 2 ), a partition Π that splits them would yield a higher Q (w,g;w ,g ) (Π) than the partition Π c that combines the two sub-populations into one leaf 1+2 (a simple proof can be found in Appendix A).
To avoid using the same information for selecting the partition and for the estimation, Athey and Imbens (2016) propose to estimate the effects in a separate sample than the one used to build the tree. They call this an 'honest' causal tree. We denote by V est the estimation set and by V tr the training set (or discovery set) that is used to build the tree (by evaluating the splitting function).
The training set and the estimation set are here obtained by taking two random subsets K tr and K est of the clusters in K. Hence, V est = k∈K est V k and V tr = k∈K tr V k = V /V est . This random split of the sample avoids any dependencies between training and estimation sub-samples. In this 'honest' version, the splitting function can be estimated as follows: where N est = |V est |. The proof follows from Athey and Imbens (2016). In (15) we can see that the splitting function is such that splits will be chosen so to maximize the heterogeneity across leaves as well as to minimize the average variance in the estimated effects. The idea is to identify the most heterogeneous partitions while introducing a penalization term which corrects the objective function to minimize the leaf-specific variation in the estimated effect. This penalization term has also the effect of reducing the depth of the tree, because leaves with a small number of observations N ( ) will exhibit a higher variance. Hence, the final depth of the tree will depend on the sample size as well as on the randomization design and the exposure mapping function. In addition to this penalization, we will also add a stopping rule based on a minimum number of observations for each condition (w, g) and (w , g ) that we are comparing. This stopping rule is required to avoid having leaves where the effect τ (w,g;w ,g ) cannot be estimated because there are no observations with observed treatment (W ik , G ik ) equal to the values (w, g) or (w , g ).

Composite splitting criterion
We now introduce a composite splitting rule targeted to multiple causal estimands. When interference is in place targeting strategies might involve both treatment and spillover effects. For example, in settings with limited resources the treatment should be provided to those who would benefit from it, i.e., with a non-zero treatment effect, whereas we could save resources by not giving the treatment to those who would benefit from other people being treated, namely units with high spillover effects. For instance, this is the case in marketing interventions where we can provide advertisements only to those who would be affected and who are less likely to get the information from someone else. Another interesting example can be found in the potential challenges of the COVID-19 vaccine distribution. Those at high risk of getting infected, even if those in close contact were immune, should be targeted. These would be individuals who are often in crowded spaces, either in the workplace or in a social setting, and would likely get infected in these environments. Therefore, these individuals are characterized by a high treatment effect even with all close contacts being treated. On the contrary, those who are in contact with a low number of people and can greatly gain from having one of these contacts vaccinated could be left without the vaccine, at least in the early stages of the distribution.
For these kind of targeting strategies involving more than one causal effect we must partition the population into sub-groups that show high level of heterogeneity in all estimands of interest.
Building a separate tree for each causal effect would provide us with different partitions that cannot be used for the design of multi-effect strategies. Therefore, we propose a composite splitting function that would result in a tree that maximises heterogeneity in all the causal estimands of interest.
This composite objective function is a weighted average of the effect-specific splitting functions: Q T (Π) = (w,g;w ,g )∈T γ (w,g;w ,g ) Q (w,g;w ,g ) (Π) with γ (w,g;w ,g ) = ω (w,g;w ,g ) τ (w,g;w ,g ) 2 where ω (w,g;w ,g ) ∈ [0, 1] is a customized weight for each estimand andτ (w,g;w ,g ) is the estimated effect in the whole sample. Each effect τ (w,g;w ,g ) , where (w, g; w , g ) ∈ T , contributes to the global objective function according to a specific weight γ (w,g;w ,g ) . γ (w,g;w ,g ) is proportional to a customized weight ω (w,g;w ,g ) , which is set by the researcher according to the extent to which the estimand τ (w,g;w ,g ) is of interest, and is normalized by the the estimated effect in the whole sample to rule out any dependence on the magnitude of the effect. The composite criterion requires that at least two of the four weights are strictly greater than zero. A similar composite objective function can be derived from the splitting functions for the 'honest' causal trees: Q H T (Π) = (w,g;w ,g )∈T γ (w,g;w ,g ) Q H (w,g;w ,g ) (Π) with γ (w,g;w ,g ) = ω (w,g;w ,g ) τ (w,g;w ,g ) 2 . (17)

Network Causal Tree Algorithm
Compared with the standard HCT algorithm, the main novelties of NCT are the introduction of interference and the possibility of including more than one effect. Specifically, the extent to which each effect τ (w,g;w ,g ) , with (w, g; w , g ) ∈ T , contributes to the determination of the tree is specified by the weight w(w, g; w , g ). Here we describe the key steps of the NCT algorithm, including the recursive partitioning based on the splitting functions and the stopping rules.

Key steps of the NCT algorithm
The proposed algorithm takes mainly six elements as inputs.
• The sample V , which collects for each unit ik the individual treatment assignment status W ik , the observed outcome Y ik and a vector of characteristics X ik ; • The network information, which is fully described by the global adjacency matrix A, including the cluster-specific blocks A k ; • The specification of the exposure mapping function g(·) which together with the adjacency matrix and possibly covariate matrix will be translated in the computation of the observed network exposure G ik for each unit; • The experimental design which will determine the computation of the probabilities π ik (w, g) and π ik (w, g, w , g ).
• The fifth aspect pertains to the weight ω w,g,w g for each causal effect; • The specification of the two parameters: maximum depth, that is the maximum depth of the tree, and the minimum size, that is the minimum number of units falling in each exposure condition (w, g) in each leaf.
After some preliminary steps, the algorithm consists of two main phases. The first phase is focused on the selection of the partition, i.e. the tree, while the second phase is concerned with the estimation of causal effects and returns point estimates and standard errors of the conditional average causal effects, for all the comparisons of interest and within each leaf of the detected partition. We report below the key steps of the NCT algorithm: 0. Phase 0 (Preliminaries): In a preliminary stage the algorithm computes the quantities and tools that will be used in the subsequent steps. In particular: (a) Given the adjacency blocks A k and potentially the covariate matrix X, for each unit the network exposure variable G ik is computed according to the rule expressed in Assumption 3.
(c) Finally, the algorithm randomly splits the clusters between the training set V est and the estimation set V est . 8 1. Phase 1 (Tree Discovery): the first step of the algorithm sprouts the tree structure of the NCT, that is, it detects the relevant heterogeneous partitions. Note that this step is performed over the discovery set only. In particular, the NCT algorithm works with the clusters belonging to the set K tr and builds the tree using a binary recursive partitioning.
(a) Recursive Partitioning. the algorithm grows a tree by maximizing the in-sample splitting criterion at each binary split. At iteration r−1 the partition can be represented as follows: ..,P is the feature that was split at iteration m and The variable x m split at iteration m together with the cutoff point c m compose a node of the tree. At iteration r, the partition will be complemented with a split of a variable x r ∈ {x p } p=1,...,P /{x m } m=1...,r−1 at some cutoff point c r : Among all the candidate splits x r and c r , the algorithm will choose the one that maximizes the in-sample splitting function in (14)

Simulation Study
Our algorithm provides an interpretable method to detect and estimate heterogeneous effects in the presence of clustered network interference. In this section we evaluate, through a set of simulations, the performance of the proposed algorithm with respect to both discovery and estimation. In particular, we investigate its ability to correctly identify the actual heterogeneous sub-populations, comparing the use of single or composite splitting functions, and we assess the performance of the Horvitz-Thompson estimator for leaf-specific treatment and spillover effects. While the latter performance assessment is quite standard in the literature, the former is critical for the development of interpretable algorithms for heterogeneous causal effects Lee et al., 2020). We evaluate the performance of the algorithm and the estimator in settings that differ Algorithm 1 Overview of the NCT algorithm • Inputs: i) observed data {W ik , Y ik , X ik } ik∈V ; ii) global adjacency matrix A, which comprises the cluster-specific blocks A k ; iii) experimental Design; iv) vector of weights ω(w, g; w , g ), where (w, g; w , g ) ∈ T ; v) tree parameters: maximum depth and minimum size.
• Outputs: (1) a partition Π if the covariate space, and (2) point estimates, standard errors and confidence intervals of the conditional average causal effects: 1. Phase 0 (Preliminaries): compute G ik and both the marginal and joint exposure probabilities π ik (w, g) and π ikjk (w, g; w , g ). Then, randomly assign clusters to discovery and estimation samples.
2. Phase 1 (Tree Discovery): using the discovery sample, build a tree according to the insample splitting criterion and stop when either the tree has reached its maximum depth or any additional split would generate leaves, which are not sufficiently representative of the four exposure conditions.
3. Phase 2 (Estimation): use the Horvitz-Thompson estimator on the estimation sample to estimate the leaf-specific CACE and their standard errors in each leaf.
with respect to three main factors: (1) the structure of the heterogeneity, (2) the extent of the effect heterogeneity, and (3) the number of clusters. In Appendix C, we also consider two additional factors: (4) the correlation structure in the covariate matrix and (5)  All simulations are performed under Bernoulli trials, that is treatment is randomly assigned independently to each unit with a fixed probability π W ik = α. In our simulation study, we also assume that interference only takes place at the neighborhood level and we choose the following definition of the network exposure: that is, the network exposure of unit ik is 1 if at least one neighbor is treated. A binary network exposure together with a binary individual treatment results in a joint treatment with four categories -i.e., (W ik , G ik ) ∈ {(w, g) = (0, 0), (1, 0), (0, 1), (1, 1)}. The binary definition of the network exposure is chosen to allow the growth of deeper trees. Given that the minimum size requirement stops the algorithm when the number of units in a child leaf is not enough to estimate a conditional causal effect, a joint treatment with four categories ensures that this stopping condition is unlikely to be met during the first few splits.
In addition, the assumption of neighborhood interference allows the computation of the marginal and joint probabilities without the need for intensive estimation procedures. In fact, the approximate algorithm for estimating the marginal and joint probabilities proposed by Aronow and Samii (2017) is computationally demanding and could not be incorporated in our simulation study. However, in the case of Bernoulli trial and network exposure defined as in (18), the probability π ik (w, g) can be computed using the formula in (3) while the joint probability π ikjk (w, g, w , g ) is simply the product of π ik (w, g) and π jk (w , g ) if the two units ik and jk are independent. On the contrary, if N wg ik and N wg jk = 0 overlap the joint treatment of the units ik and jk will be dependent, that is, π ikjk (w, g; w , g ) = π ik (w, g)π jk (w , g ). In this case, the joint probability can still be readily computed using combinatorics formulas on two overlapping sets.

Data generating process
For each simulation m = 1 . . . , M we generated K clusters and simulated, within each cluster, Erdős-Rényi random graphs (Erdős and Rényi, 1959) with n k = 100 nodes and a fixed probability (0.01) to observe a link. Given the definition of the network exposure we removed isolated nodes from the analysis to make the Assumption 5 hold. W ik and any of the 10 covariates X ip were sampled from independent Bernoulli distributions with probability 0.5: W i ∼ Ber(0.5) and X ip ∼ Ber(0.5).
In the simulation study we focus on two main effects: (i) the pure treatment effect τ (1,0,0,0) , and (ii) the pure spillover effect τ (0,1,0,0) . To ease notation, we denote by τ the treatment effect and by δ the spillover effect. After setting the value of these two conditional effects for each unit with covariates X ik (depending on the simulation scenarios), we generated the four different potential outcomes as follows: Finally, the observed outcome is given by: We now detail how we varied the three factors (1), (2), and (3). We simulated two different scenarios with respect to the heterogeneity structure (1). In the first scenario we have: Hence, in this scenario the heterogeneity driving variables (HDV), i.e., X i1 and X i2 , are the same for both the treatment effect τ and the spillover effect δ and the two causal rules overlap. In the second scenario, we introduce a change in the drivers of the heterogeneity in the following way: Hence, in the second scenario the heterogeneity drivers are different for the two causal effects.
Specifically, we have: X i1 and X i2 for the treatment effects, and X i1 and X i3 for the spillover effects. In addition, we have two causal rules for the treatment effect, namely {X ik1 = 0, X ik2 = 0} and {X ik1 = 0, X ik2 = 1} and two different causal rules for the spillover effect, namely {X ik1 = 1, X ik3 = 0} and {X ik1 = 1, X ik3 = 1}. For each structural scenario we varied the effect size: h = (k + 0.1) 10 k=1 . Figure 4 graphically represents the two simulations' scenarios.

Yes No
Yes No (a) First scenario.

Performance measures
In all the scenarios the performance of the NCT algorithm is evaluated using the following measures, averaged over the M generated data sets: • Average number of correctly discovered heterogeneous causal rules corresponding to the leaves of the generated NCTs (reported with respect to the effect sizes in Figures 5 and 6) in the discovery sample; • Average conditional treatment and spillover effects estimated for each correctly detected heterogeneous terminal leaf (reported asτ andδ in Tables 1, 2, 3 and 4); • Average standard errors estimated for each correctly detected heterogeneous terminal leaf (reported asŝe(τ ) andŝe(δ) in Tables 1, 2, 3 and 4); • Monte Carlo bias in the estimation sample: where Π m is the partition selected in simulation m; • Monte Carlo MSE in the estimation sample: • Coverage, computed as the average proportion of units for whom where the estimated 95% confidence interval of the causal effect in the assigned leaf includes the true value:

Results
We will first analyze the ability of the algorithm to correctly detect the heterogeneous subgroups in the first simulation scenario, that is, when the heterogeneity is the same for the two causal effects of interest. Figure 5 reports the average number of correctly discovered heterogeneous causal rules with composite splitting rule or effect-specific splitting rules targeted to the treatment effect or the spillover effect, in the case of 10, 20 and 30 clusters. Given that in this scenario the causal rules defining the heterogeneity is the same for both effects, all the splitting rules, targeted to a single effect or to both effects, have a similar performance and are able to detect the two heterogeneous leaves with a success rate that gets higher as the effect size increases. In addition, as the number of cluster grows the minimum effect size allowing the algorithm to optimally discover all the heterogeneous sub-populations gets lower.
Tables 1, 2, 3 report the results for the first scenario 4 for the performance of the estimator in the correctly detected leaves. We only report the results of the estimation procedure on the tree built with the composite splitting rule, as the spitting rule would only affect the identification of the heterogeneous sub-populations but not the estimation of the causal effects once these subpopulations are correctly detected. The estimator is able to estimate the heterogeneous treatment and spillover effects without bias and its precision grows as the number of clusters increases. Interestingly, NCT provides more accurate estimates of the heterogeneous spillover effects than the treatment effects. This is due to the larger number of units with (W ik = 0, G ik = 1) than those with , by definition of the network exposure (by design, units with (W ik = 0, G ik = 1) are on average 50% more than units with (W ik = 1, G ik = 0)). The higher precision in the estimation of the spillover effects is also reflected in the identification of the heterogeneous subgroups, which is more accurate when splits are targeted to the minimization of the MSE of the spillover effect ( Figure 5).

Empirical Application
In this section we provide an empirical application of the proposed methodology. We use data from a randomized experiment designed to assess the effectiveness of intensive information sessions to promote the uptake of a new weather insurance policy among farmers in rural China. The promoted policy is addressed to rice farmers and it is aimed at protecting their products from adverse weather shocks (Cai et al., 2015). The sample consists of 4,569 households belonging to 47 different villages.
Households were randomly assigned to the the intensive training session. 10 Here, the individual treatment variable W ik ∈ {0, 1} for each household i living in village k represents the assignment to the intensive (W ik = 1 ) or simple (W ik = 0) information session.
Households are linked according to an observed village-specific friendship network, represented in Figure 7. Relationships between households belonging to different villages are negligible and were removed. Friendship relationships between households in a given rural village k are fully described by the adjacency matrix A k , where the generic element A k (i, j) equals 1 if the household i in village k has nominated the household j in the same village k as friend. Therefore, the adjacency matrix is not symmetric. We denote with N ik the set of nominated friends of unit ik, and with N ik the out-degree. Figure 8a shows the overall degree distribution, while Figure 8b represents the distribution of the number of treated neighbors. In this setting interference could take place because the information received in the intensive information session could be transferred to network contacts. We assume here an exposure to the network treatments only at the neighborhood level and, in particular, we define the network exposure variable as in (18). Hence, G ik is equal to 1 if a unit has at least one treated friend, and 0 otherwise. This definition of the network exposure is justifiable by the fact that the decision of a farmer to buy the weather insurance might depend only on the information either received directly in the intensive session or indirectly by one of his friends. This is also what has been found in Cai et al. (2015) where the only significant spillover effect is from first neighbors.  Table 5.  To evaluate the heterogeneity of treatment and spillover effects, we included in the analysis all the observed characteristics that could plausibly moderate the effect of the intervention: three dummies representing the household's production area (reg1, reg2 and reg3, respectively); a binary variable which equals 1 if the head of the household is at least 50 years old (age); a binary variable being 1 if the household's head has successfully completed high school (educ); a binary variable distinguishing families who are strongly worried about weather phenomena (prob dis, which is equal to 1 if perceived probability of a weather disaster happening in the coming year exceeds 0.30); a binary variable which identifies families who declare to be risk averse (averse).

Application's Results
In this section, we present the most relevant empirical findings. The following figures report the trees built targeting single effects ( Figure 10) and multiple effects ( Figure 11). In each node, we also report the estimates of the main treatment and spillover effects and the size of the selected sub-populations. 20 villages were randomly assigned to the training set, while the remaining villages have been allocated to the estimation set. 11 Furthermore, we have set the maximum depth at 3 to maintain a high level of interpretability, while we have set at 20 the minimum number of units that must be present in the child leaves for each of the four exposure conditions.
We can see that, in the whole population, the treatment has a positive effect on insurance take up, while the spillover effect is negligible. The most relevant heterogeneity drivers result to be the perceived probability of disaster and the production area (specifically, living in the third area).
However, the estimated tree slightly changes under different specifications of the spitting rule. Figures 10a and 10b separately represent the trees targeted to τ (1, 0, 0, 0) and τ (0, 1, 0, 0). The most important variable driving the heterogeneity of the treatment effect is the perceived probability of disaster. The treatment appears to be particularly effective in the sub-population including less concerned and older households, with both a high or low level of education. We can speculate that younger households who are more worried about possible disasters benefit less from an intensive information session because even less information would prompt them to purchase a weather insurance. When looking at spillover effects, living in the third region is the main characteristic determining the heterogeneity of spillover effect. In this region younger households with a lower concern about a possible disaster are those who benefit the most from receiving the information about the weather insurance from some of their friends.
11 The random assignment of the villages to the training and estimation samples has been kept fixed in all the trees.  Figure 11b is related to the tree which has been built assigning an equal weight to all the four effects in T . In this application, the composite tree which considers both τ (1, 0, 0, 0) and τ (0, 1, 0, 0) coincides with the tree based on τ (1, 0, 0, 0) only. Therefore, the latter leads the composite variability across partitions. Finally, the network causal tree that incorporates all the four effects shows slightly different results: here the sub-population where the treatment results to be more effective is the one including older households who are currently settled in the second production area.

Conclusions
Depending on our characteristics we might respond differently to a treatment or an intervention.
Similarly, we might be more or less susceptible to the influence of other people who have received the treatment. Understanding the heterogeneity of the effect of a treatment with the aim of targeting people who would benefit from it has been the focus of a recent field of research, especially applied to medicine. Investigating how different people respond differently to the treatment received by others can be crucial, particularly in settings with limited resources where spillover effects could be leveraged. In this paper, we have introduced a new algorithm to estimate heterogeneous causal treatment and spillover effects in the presence of clustered network interference. The proposed network causal tree model bridges the gap between two streams of causal inference literature: estimators for causal effects under interference and tree-based methods for the discovery of heterogeneous sub-populations. We build upon the seminal algorithm proposed by Athey and Imbens (2016) to account for clustered network interference through a rework of the criterion function. Leaf-specific causal effects are then estimated using the Horvitz-Thompson estimator proposed by Aronow and Samii (2017).
The proposed NCT algorithm has enhanced interpretability and shows an excellent performance in a set of Monte Carlo simulations. In particular, the algorithm is able both to spot the relevant sources of heterogeneity in the data and to consistently estimate the conditional treatment and spillover effects. Moreover, we have introduced a composite splitting function that allows the researchers to simultaneously detect the sub-populations where both treatment and spillover effects are heterogeneous. The identification of these multi-effect heterogeneous subgroups is crucial for the design of targeting strategies that involve multiple effects. For instance, in marketing campaigns a person might not be affected by an advertisement received directly by a company but might be susceptible to the advertisement received by her friends. In this case, resources could be saved by promoting the product among people who could be directly susceptible and letting those who could be more influenced by their friends receive the advertisement indirectly. Our simulation study shows that the use of such composite splitting rule is able to correctly detect all the heterogeneous sub-populations defined by both treatment and spillover effects, as well the ones defined by one effect only. Therefore, the selected partition of the population can be used to design strategies whose objective function incorporates multiple effects, but can also be used a posteriori to target subgroups maximizing either a treatment or a spillover effect.
When applied to real-world data, the NCT algorithm provided useful insights on the effectiveness of intensive training sessions among Chinese rural households on the uptake of a weather insurance policy. We found that the main characteristics responsible for the heterogeneity of the effects are perception of a possible disaster, the production area, and the age of the farmers. However, these heterogeneity drivers play a different role with respect to the two main treatment and spillover effects and when the tree is built using composite splitting function.
Nevertheless, the proposed algorithm may suffer from the limitations common to tree-based methods: instability to the random allocation of units in the training sample and the potential impact of outlier observations in the node-specific estimations. Here, we propose an algorithm that selects a single tree because of its high interpretability that plays a fundamental role in policy design. Indeed, the single tree algorithms are suitable for the discovery of heterogeneous subpopulations, which can give useful insights into the main variables that drive the heterogeneity and can be used to design targeting strategies. Moreover, we did not find instability in neither the detection nor the estimation of heterogeneous effects in the Monte Carlo simulations. Nonetheless, the proposed algorithm could be extended to tree-based ensemble methods, following Wager and Athey (2018) and . By averaging the estimates from many single trees, this extension could enhance estimation precision at the cost of reduced interpretability (see Lee et al. (2020) for a discussion on the trade-off between accuracy and interpretability). In addition, our approach might be rearranged to deal with settings where the network cannot be partitioned into well-defined and pre-specified clusters. However, in some network structures, clusters could be detected by implementing a network-based community detection algorithm (Fortunato, 2010) and the NCT algorithm could be applied on the detected communities, while the estimator should take into account the uncertainty in group membership .
Here we assume that the network exposure variable is discrete, with the performance of the algorithm being affected by the number of categories resulting from the exposure mapping function.
In our simulation study and application we have used a binary neighborhood exposure, which allowed us to grow deeper trees, reduce the number of possible causal effects, and have enough observations for each exposure condition to maintain the variance in a reasonable range. However, alternative specifications could be used. For instance, the network exposure could be defined as the proportion of treated neighbors, perhaps categorized into few bins. A more complex definition of the network exposure, possibly resulting in a continuous variable, would require some methodological adjustments in the estimation strategy, but it would allow to model a wider ensemble of real-world interference mechanisms. We leave this extension to future work.
Finally, further research is needed to use the selected partition to actually design targeting strategies involving both treatment and spillover effects. Furthermore, these strategies should also rely on the average susceptibility of network contacts as well as on heterogeneous influential power.
To this end, new causal estimands could be defined and their heterogeneity should be investigated and incorporated in the design of complex targeting rules, aimed at maximizing the effect of the intervention while saving resources.
where the expectation is over the randomization distribution of W ik and the induced distribution on G ik and the first equality holds by consistency.
Proposition 6 (Population Unbiasedness). The estimator is unbiased with respect to the population mean of the potential outcomes: where the expected value is taken over the sampling distribution. Proof.
The proof of unbiasedness of τ (w,g,w ,g ) follows directly.
Proof. As in Proposition 1 µ (w,g)) ( (x)) is unbiased. Hence, for consistency to hold we need to prove that the variance goes to 0 as N goes to infinity. Following Aronow and Samii (2017), it is easy to show that the variance of µ (w,g) ( (x)) is given by: Since |Y ik (w, g)|/π ik (w, g) ≤ C < 1 and given that in each cluster the sample size belonging to leaf (x) is bounded, i.e., n k ( (x)) ≤ B( (x)) ≤ B, we have: Hence, Q (w,g;w ,g ) (Π) ≥ Q (w,g;w ,g ) (Π c ).

C Additional Monte Carlo Simulations
We included in the simulation study two additional sets of simulations where (4) we introduce correlation between the covariates, and (5) we replace the Erdős-Rényi model for network formation with an exponential random graph (ERGM) model introducing homophily within the clusters. These two additional simulations are conducted with 30 clusters and under the first scenario introduced in Section 5, where the heterogeneity is the same for the two causal effects of interest.

C.1 Correlated covariates
In Figure 12 we report the number of correctly detected leaves under low and high correlation (0.25 and 0.5), while in Tables 6 and 7 we report the estimated treatment and spillover effects with their standard error in the two heterogeneous leaves, together with the MSE, bias and coverage of the average treatment and spillover effects in the sample. From Figure 12 one can see that the correlation between covariates compromises the ability of the algorithm to correctly identify the heterogeneous subgroups. This is due to the fact that, as the covariates become more similar to each other, it becomes harder for the algorithm to detect the true HDV. Such a problem is common to all tree-based algorithms. Hence, we argue that one should carefully check the correlation patterns between the variables to get a sense of the reliability of the discovered subgroups.  NCT (treatment) NCT (spillover) Figure 12: Simulations' results for correctly discovered leaves in the first scenario with correlated covariates  Nevertheless, for both correlation levels (0.25 and 0.50) the estimator seems to perform well within correctly detected leaves (see Tables 6 and 7). Table 8 shows the results in the case of network homophily within the clusters. In this case, we find larger standard errors than the original scenario reported in 3 without homophily. As a consequence, the Monte Carlo MSE is also slightly larger.

C.2 Network homophily within the clusters
Moreover, as we can see from Figure 13 there is a decrease in the ability of the algorithm to discover the true leaves. Indeed, if one compares this Figure with the right panel of Figure 5, one can see how the correct discovery of the true leaves is slower in the case with homophily network.
This is due to the fact that the standard errors are larger than in the original first scenario.