# A Bayesian Approach to Estimate Causal Peer Influence Accounting for Latent
Network Homophily

Seungha Um
Division of Biostatistics
Department of Population Health
New York University Grossman School of Medicine
New York, NY 10016
<EMAIL_ADDRESS>
&Tracy Sweet
Department of Human Development and
Quantitative Methodology
University of Maryland
College Park, MD 20742
<EMAIL_ADDRESS>
&Samrachana Adhikari
Division of Biostatistics
Department of Population Health
New York University Grossman School of Medicine
New York, NY 10016
<EMAIL_ADDRESS>

###### Abstract

Researchers have focused on understanding how individual’s behavior is
influenced by the behaviors of their peers in observational studies of social
networks. Identifying and estimating causal peer influence, however, is
challenging due to confounding by homophily, where people tend to connect with
those who share similar characteristics with them. Moreover, since all the
attributes driving homophily are generally not always observed and act as
unobserved confounders, identifying and estimating causal peer influence
becomes infeasible using standard causal identification assumptions. In this
paper, we address this challenge by leveraging latent locations inferred from
the network itself to disentangle homophily from causal peer influence, and we
extend this approach to multiple networks by adopting a Bayesian hierarchical
modeling framework. To accommodate the nonlinear dependency of peer influence
on individual behavior, we employ a Bayesian nonparametric method,
specifically Bayesian Additive Regression Trees (BART), and we propose a
Bayesian framework that accounts for the uncertainty in inferring latent
locations. We assess the operating characteristics of the estimator via
simulation study. Finally, we apply our method to estimate causal peer
influence in advice-seeking networks of teachers in secondary schools and
assess whether teachers’ beliefs about mathematics education is influenced by
the beliefs of their peers from whom they receive advice. Our results suggest
that overlooking latent homophily can lead to either underestimation or
overestimation of causal peer influence, accompanied by considerable
estimation uncertainty.

_Keywords_ Causal Peer influence Effect $\cdot$ Latent Homophily $\cdot$
Causal Inference $\cdot$ Bayesian additive regression trees $\cdot$ Social
network

##  1 Introduction

Social network data typically consist of a set of individuals and the
relational ties among them. Within a social network, an individual’s behavior
may adapt or respond to the behaviors of others. This phenomenon of alters’
behaviors impacting an ego’s behavior is commonly referred to as peer
influence. Identifying and quantifying peer influence has been a focus of
research in a variety of disciplines, including substance use (Napier et al.,
1984; Henneberger et al., 2021), health outcomes (Christakis and Fowler, 2007,
2008), and education (Supovitz et al., 2010; Rosenqvist, 2018). As an example
from education research, consider an advice-seeking network among teachers in
a school. In such a network, a teacher’s individual beliefs about instruction
could be influenced by the beliefs of their advisors and other teachers from
whom they seek advice (Supovitz et al., 2010; Sweet and Adhikari, 2020).
Understanding peer influence in the network of teachers may be important for
enhancing instructional methods and facilitating the dissemination of
innovative teaching approaches.

Despite the widespread interest in peer influence, identifying and accurately
estimating a causal peer influence effect using observational network data is
often challenging due to confounding by multiple measured and unmeasured
factors that lead to homophily (Ogburn, 2018). Homophily is the tendency of
individuals to associate and form connections with others who are similar to
themselves (McPherson et al., 2001; Ogburn, 2018; Shalizi and Thomas, 2011).
In our example with the advice-seeking network, two teachers in a school may
seek advice from each other because they share similar beliefs, for example,
as a consequence of both teaching the same grade. If they both also encourage
their students to construct their own knowledge through exploration and
experimentation, reflecting similar teaching beliefs, it is important to
separate the influence of peer teaching beliefs from similarity in beliefs to
establish causality.

Thus, to accurately identify causal peer influence, the effect of peers’
behaviors should be distinguished from the effect of similarities, which may
induce similar behaviors among peers (Shalizi and Thomas, 2011). However, for
many social networks it is not possible to observe or measure all homophilous
attributes that lead to relational ties, resulting in latent homophily.
Further, when there is confounding due to latent homophily, the peer influence
effect is not causally identified without making additional assumptions
(Shalizi and Thomas, 2011). In fact, none of the previous work that focused on
estimating peer influence in advice-seeking networks for teachers addressed
confounding due to latent homophily (Spillane et al., 2018; Sweet and
Adhikari, 2020).

In this paper, we propose an approach to estimate causal peer influence
effects while adjusting for confounding due to latent homophily. Specifically,
we present identification assumptions for the causal peer influence effect in
the presence of latent homophily, extending previous work (McFowland III and
Shalizi, 2021; Cristali and Veitch, 2022), and develop a flexible Bayesian
hierarchical framework for estimation. We further extend our proposed Bayesian
hierarchical framework to estimate peer influence effect from multiple
networks while accommodating the nonlinear dependency of peer influence on
individual behavior.

###  1.1 Related Work

A common strategy to account for unobserved confounders in causal inference
involves the use of proxies of unmeasured confounders (Liu et al., 2024).
Proxy variables, commonly used in epidemiology studies and referred to as
negative control variables, are employed to detect and sometimes correct for
unmeasured confounding (Tchetgen, 2014; Miao et al., 2018; Shi et al., 2020).
There are two primary approaches to leverage proxy variables: i) selecting a
proxy variable from a set of plausible observed candidates, or ii) estimating
a proxy variable as a latent variable using observed data. This strategy of
using proxy variables is applicable regardless of whether individuals are
independent or connected through network ties, with further elaboration
provided below.

The first approach of selecting a proxy variable leverages the association
between an unmeasured confounder and a proxy variable, enabling researchers to
establish the conditional distribution of the proxy variable given the
unmeasured confounder. Recent papers employing such a proxy variable have
focused on the identification and estimation of causal effects involving both
non-network data (Liu et al., 2024; Miao et al., 2018; Wang Miao and Zhou,
2023; Louizos et al., 2017) and network data (Egami and Tchetgen Tchetgen,
2023; Liu and Tchetgen Tchetgen, 2022). For example, Louizos et al. (2017)
viewed the observed proxy variables as noisy versions of unmeasured
confounders and proposed a method that was robust against varying noise
levels. In considering network data to account for latent homophily, Egami and
Tchetgen Tchetgen (2023) used two proxy variables; one was known not to be
causally affected by the treatment of interest, and the other did not causally
affect the outcome of interest.

However, identifying suitable proxy variables from observed candidates in
real-world scenarios can often be challenging. To address this challenge, the
second approach of estimating proxy variables may be more appropriate. In a
network context, we can leverage information embedded within network ties to
estimate an unmeasured homophilous attribute. In fact, many existing latent
variable social network models (Dabbs et al., 2020) can be utilized to
estimate such proxy variables. For example, consider a latent space network
model (Hoff et al., 2002) which assumes that each node in a network can be
represented by a position in a lower dimensional latent space. Conditional on
the latent positions, the probability of two nodes sharing a tie is inversely
related to the distance between their latent positions in this space.
Consequently, these latent positions, which capture the existing network
structures including homophily, transitivity and clustering, can serve as a
proxy for latent homophilous attributes. Adjusting for these inferred latent
positions, McFowland III and Shalizi (2021) introduced a consistent estimator
of causal peer effects. Alternatively, Cristali and Veitch (2022) used
embedding methods to capture the latent properties of each node as proxies and
provided a nonparametric identification of causal peer influence. However,
both studies are limited to investigations assessing a linear relationship
between ego’s and peers’ behaviors, where a coefficient in a linear model was
interpreted as the causal peer influence effect (McFowland III and Shalizi,
2021; Cristali and Veitch, 2022). Such restrictive assumptions of linearity
presents challenges for modeling nonlinear relationships between behavior and
peer influence, including interactions between predictors. Also, both of these
studies on the estimation of causal peer influence are examined under a single
network setting. Therefore, a comprehensive and flexible causal framework for
assessing peer influence in multiple social networks is lacking.

###  1.2 Motivating application

ID | size | degree | density | transitivity | score (2012) | peers’ score (2012) | score (2013) | difference
---|---|---|---|---|---|---|---|---
S1 | 27 | 3.1(1.4) | 0.2 | 0.34 | 19.8(2.6) | 19.8(1.3) | 19.9(2.5) | 0.9(0.7)
S2 | 27 | 3.6(3.0) | 0.3 | 0.32 | 18.4(2.2) | 18.5(3.8) | 19.4(1.6) | 1.8(1.5)
S3 | 25 | 3.0(1.1) | 0.3 | 0.42 | 19.0(2.4) | 19.1(1.4) | 19.5(2.6) | 1.6(1.6)
S4 | 30 | 3.5(1.3) | 0.2 | 0.51 | 19.5(2.1) | 19.1(3.8) | 20.2(2.3) | 1.5(1.0)
S5 | 23 | 2.2(1.3) | 0.2 | 0.31 | 19.3(2.4) | 18.9(4.6) | 19.2(1.8) | 1.8(1.5)
S6 | 26 | 2.3(1.0) | 0.2 | 0.46 | 17.8(2.3) | 17.1(3.8) | 18.6(2.2) | 1.7(1.4)
S7 | 20 | 2.3(1.7) | 0.2 | 0.37 | 17.7(2.3) | 16.0(5.7) | 19.0(1.7) | 2.0(1.7)
S8 | 26 | 3.4(1.5) | 0.3 | 0.45 | 18.8(2.1) | 18.7(1.2) | 19.7(1.7) | 1.5(1.3)
S9 | 24 | 2.6(1.1) | 0.2 | 0.14 | 18.2(2.7) | 18.1(2.1) | 16.9(1.4) | 2.2(2.2)
S10 | 28 | 2.8(1.3) | 0.2 | 0.39 | 19.1(1.8) | 18.3(4.0) | 18.0(1.6) | 1.9(1.9)
S11 | 23 | 2.2(1.1) | 0.2 | 0.21 | 19.1(2.6) | 17.8(4.4) | 18.9(2.2) | 1.8(1.8)
S12 | 22 | 2.2(0.7) | 0.2 | 0.55 | 18.7(2.3) | 18.8(1.3) | 19.6(2.2) | 2.3(2.4)
S13 | 15 | 2.6(1.9) | 0.4 | 0.18 | 18.8(2.4) | 18.1(5.3) | 19.1(2.1) | 1.9(1.3)
S14 | 13 | 2.5(1.6) | 0.4 | 0.16 | 19.8(2.0) | 17.3(7.8) | 19.8(2.4) | 1.4(1.3)
Table 1: Summary statistics of the network characteristics and teaching belief
scores for each school with mean (sd). The difference represents the absolute
difference between the scores from 2012 and 2013.

Returning to our example of advice-seeking networks, our work is motivated by
a real-world scenario involving elementary school staff in a mid-sized
Midwestern school district referred to as Auburn Park (Pitts and Spillane,
2009). In 2012, this suburban district served approximately 5,800 students
across 14 elementary schools. In Spring 2012 and Spring 2013, staff members in
each school responded to surveys that included items related to their beliefs
about mathematics education and about to whom they sought advice and
information about curriculum and instruction. In their study on mathematics
beliefs, Spillane et al. (2018) identified a mathematics-related belief
construct that they called student-centered learning; this was based on 5
Likert scale items that advocated for engaging students with guiding questions
rather than direct instruction, promoting an environment where students can
construct their own knowledge through exploration. They used the sum score to
generate an overall belief score for each respondent, which we will use as the
outcome of interest in our study. Details on the specific survey items can be
found in the Appendix (Table 2) and summary statistics of the network
characteristics and teaching belief scores for each school are represented in
Table 1. In Table 1, for each school we present the size of the network
(size), average degree (average number of ties per node), network density (the
ratio between the actual number of ties and the maximum possible number of
ties), transitivity (the ratio between the observed number of closed triplets
and the maximum possible number of closed triplets), along with a summary of
the belief scores in 2012 and 2013. The advice-seeking networks observed
across the 14 schools varied in their size, ranging from 13 to 30 teachers,
and also in their network structures with network densities ranging from 0.18
to 0.41. Further information about the dataset is available in Spillane et al.
(2018) and Sweet and Adhikari (2020).

We can visualize changes in teachers’ individual beliefs from 2012 to 2013,
potentially in response to their peers’ belief scores in 2012. As an
illustration, Figure 1 shows the belief scores for advice-seeking networks in
School 8 and School 9; similar plots for the remaining schools are shown in
Figure 11 in the Appendix. For effective visualization, we dichotomized the
continuous math belief score into high and low categories using a cut-off of
$18.85$, the average belief score from 2012. Red nodes are nodes with low
belief scores, and blue nodes are those with high belief scores.

As shown in Figure 1 (left), 17 teachers (65.3%) in School 8 exhibit high
belief scores in 2012, and 23 teachers (88.4%) have high scores in 2013. This
suggests a plausible shift towards higher scores in 2013. Conversely, for
School 9 (right), 15 teachers (57.5%) had low scores in 2012, while 22
teachers (84.6%) had low scores perhaps in 2013. These observations suggest
that individuals tended to have lower scores in 2013 when their peers had low
scores in 2012, and higher scores when their peers had high scores, which
suggests some evidence of peer influence.

Our primary goal in this paper is to assess whether the teachers in each
school changed their beliefs about mathematics education based on the beliefs
of their peers from whom they receive advice. In order to establish that the
belief scores of the teachers in 2013 were causally influenced by their peers
scores in 2012, we need to rule out all other common causes. However, this
assumption is likely to be violated in advice-seeking networks because
teachers often seek advice from peers who share similar attributes and
beliefs, which can also contribute to changes in their belief scores. Many of
these attributes are, however, not fully observed. Estimates of causal peer
influence effects that do not adjust for these unobserved homophilous
attributes are likely going to be biased. Therefore, to isolate the influence
of peers from the effects of these unobserved attributes, we propose a
Bayesian approach to estimate and adjust the latent homophilous variable.

Figure 1: The teachers’ scores about student-centered beliefs in 2012 and 2013
for school 8 and school 9 with network ties. The belief scores are represented
visually by coloring each node based on high/low belief scores, using a cut-
off of $18.85$ for effective visualization.

###  1.3 Contribution

This paper presents a novel Bayesian framework for estimating causal peer
influence while incorporating latent variables to address unmeasured
confounding due to homophily with a semiparametric approach. The framework is
designed for scenarios where independent yet similar networks are observed
simultaneously, and there is a nonlinear relationship between individual
behavior and peers’ behaviors. Our main contributions are as follows.

First, we address the issue of disentangling homophily from peer influence by
leveraging latent locations inferred from the network itself, instead of using
observed proxy variables. This extends the work of McFowland III and Shalizi
(2021) to integrate the estimation of peer influence within a causal framework
by clearly outlining causal identification assumptions with latent homophily.
Furthermore, our work focuses on inferring latent locations from the
continuous latent space model which can be viewed as a generalized version of
the latent community model considered in simulation studies by McFowland III
and Shalizi (2021) .

Second, in contrast to the existing literature, which has largely focused on
estimating causal peer influence within a single network (McFowland III and
Shalizi, 2021; Cristali and Veitch, 2022; Egami and Tchetgen Tchetgen, 2023),
we extend the latent locations strategy to multiple networks through a
Bayesian hierarchical modeling framework (Sweet and Adhikari, 2020). This
framework facilitates the estimation of shared homophilous attributes,
represented as latent locations, across different networks. Additionally, it
allows us to quantify the uncertainty in inferring these latent positions,
which is particularly useful when each observed network has a small number of
nodes.

Lastly, we adopt BART, a Bayesian nonparametric model, to address the
limitations of linear models in capturing complex, non-linear relationships
between individual behavior and peers’ behaviors while also facilitating
uncertainty quantification using posterior predictive distributions. Existing
research has primarily focused on estimating causal inference using linear
regression models, often overlooking the aspect of uncertainty quantification
(McFowland III and Shalizi, 2021; Egami and Tchetgen Tchetgen, 2023;
Forastiere et al., 2021). Although Forastiere et al. (2022) employed penalized
splines to accommodate nonlinear dependencies between outcomes and
neighborhood treatments without homophily, BART efficiently captures
nonlinearity and offers the additional flexibility of not requiring a
predefined functional form. To the best of our knowledge, this is the first
attempt to incorporate the estimation of homophily attributes within a
Bayesian framework for identifying causal peer influence and extend it to the
multiple networks setting.

The remainder of this article is organized as follows. In Section 2, we
introduce notation and establish the nonparametric identification of causal
peer influence, outlining the necessary conditions. In Section 3, we propose
estimation strategies for the causal peer influence effect within a Bayesian
framework. We then evaluate our methods through a simulated study in Section
4, and we apply our proposed methodology to the motivating advice-seeking
network dataset in Section 5. Section 6 concludes with a discussion.

Baseline Follow-up
$X_{i}$$U_{i}$$U_{j}$$X_{j}$$Y_{i1}$$Y_{i2}$$Y_{j1}$$Y_{j2}$$A_{ij}$ Figure 2:
The graphical causal model for our setting. The shaded (unshaded) ones
represent observed (unobserved) variables. The solid lines indicate causal
relations, while dotted lines indicate the peer effect of interest. For
simplicity, we illustrate causal relationship between two connected nodes.

##  2 Causal Inference in the Presence of Latent Homophily

###  2.1 Notation

We denote a network as $\mathcal{A}_{n}=(V_{n},A)$, where $V_{n}$ is a set of
nodes with cardinality $n$ and $A$ is a $n\times n$ adjacency matrix such that
$A_{ij}=1$ represents the presence of a link from node $i$ to node $j$ for
$i\neq j\in V_{n}$, and $A_{ij}=0$ otherwise. A network is undirected if
$A^{T}=A$; otherwise the network is directed. We use $\mathcal{N}_{i}$ to
denote a set of neighbors that share a tie with node $i$. We assume that a
network is fixed over time and the behavior of node $i$ is observed at
baseline $(t=1)$ and a single follow-up $(t=2)$. The behavior of node
$i\in\\{1,\cdots,n\\}$ at time $t\in\\{1,2\\}$ is denoted by $Y_{it}$ for node
$i$. To measure the peer influence on the behavior of node $i$ at follow-up
($Y_{i2}$), it is possible to consider a multivariate variable of peers’
behaviors, $\\{\bm{Y}_{j1}:j\in\mathcal{N}_{i}\\}$ for node $i$. However, in
practice, the peer influence is assumed to be explained through a specific
function $\phi(\cdot)$ designed for dimension reduction, known as an exposure
mapping (Aronow and Samii, 2017). We define the summary of peers’ behaviors at
baseline for node $i$ as
$G_{i}=\phi(\bm{Y}_{j1}:j\in\mathcal{N}_{i})\in\mathbb{R}$. A common choice
for $G_{i}$ is the average of peers’ behaviors at $t=1$, that is,
$G_{i}=\sum_{j\in\mathcal{N}_{i}}Y_{j1}/\sum_{j}A_{ij}$ (Forastiere et al.,
2021; Egami and Tchetgen Tchetgen, 2023). Alternatively, we can consider the
weighted mean of peers’ behaviors as
$G_{i}=\sum_{j\in\mathcal{N}_{i}}w_{ij}A_{j1}/\sum_{j\in\mathcal{N}_{i}}w_{ij}$,
where each peer has a distinct level of impact on the behavior of node $i$.
Throughout the paper, we consider the average of peers’ behaviors to define
$G_{i}$ for node $i$.

Let $\operatorname{\bm{X}}_{i}$ denote a vector of observed covariates for
node $i$, and $U_{i}$ denote an unmeasured confounder that affects not only
behaviors at baseline ($Y_{i1}$) and follow-up ($Y_{i2}$), but also the
network ties $A_{ij}$ for $j\in\mathcal{N}_{i}$. The corresponding directed
acyclic graph (DAG) is represented in Figure 2. To estimate the causal peer
influence of node $j$ on node $i$ using observational data, we condition on
the network tie $A_{ij}=1$ (i.e. presence of network ties). However, when
conditioning on the collider $A_{ij}=1$ , there is an unblocked backdoor path
$Y_{j1}\leftarrow U_{j}\rightarrow A_{ij}\leftarrow U_{i}\rightarrow Y_{i2}$
where $U_{j}$ and $U_{i}$ are unmeasured (Shalizi and Thomas, 2011; Egami and
Tchetgen Tchetgen, 2023). Therefore, in the presence of unmeasured confounder
$U_{i}$, identification and estimation of the causal peer influence is not
feasible without additional assumptions. We outline the identification
assumptions that rely on the latent homophily property of a network in Section
2.3. More specifically, we leverage on the observation that when a network is
formed based on homophily, it is likely that a node will be similar to its
neighbors, indicating that the values of $U_{i}$ and $U_{j}$ are close.

###  2.2 Causal Estimand: Causal Peer Influence

Using the potential outcomes framework (Splawa-Neyman et al., 1990; Rubin,
1974), we define $Y_{i2}(G_{i}=g)$ as the potential outcome for node $i$’s
behavior at follow-up. The potential outcome, simplified as $Y_{i2}(g)$,
denotes the behavior of node $i$ that would have been observed at follow-up if
the average of peers’ behaviors at baseline were $g$. To measure the influence
of peer behavior on individual’s behavior, we can define the average causal
peer influence (ACPI) given the observed network $\mathcal{A}_{n}$, denoted by
$\tau(g,g^{\prime};\mathcal{A}_{n})$, using the potential outcome notation as

$\displaystyle\tau(g,g^{\prime};\mathcal{A}_{n})=E\left[Y_{i2}(g)-Y_{i2}(g^{\prime})|i\in\mathcal{A}_{n}\right]=\frac{1}{n}\sum_{i\in\mathcal{A}_{n}}\left[Y_{i2}(g)-Y_{i2}(g^{\prime})\right]=\mu(g)-\mu(g^{\prime})$

for $g,g^{\prime}\in\mathcal{G}$ where $\mathcal{G}$ is the support of
$G_{i}$. The causal estimand is written as the empirical mean (Forastiere et
al., 2021, 2022; Egami and Tchetgen Tchetgen, 2023), which can be viewed from
a "super-population" perspective in terms of expectation (Imbens and Rubin,
2015; Forastiere et al., 2021). The expectation is conditional on the observed
network $\mathcal{A}_{n}$ and we write
$\tau(g,g^{\prime};\mathcal{A}_{n})=\tau(g,g^{\prime})$ omitting dependency on
$\mathcal{A}_{n}$ for notational simplicity.

This estimand compares the average potential outcomes at two distinct
hypothetical peers’ behaviors $g$ and $g^{\prime}$, while preserving the
identical set of individuals connected by the same network ties. The ACPI
depends on the quantity of the average potential outcome at level $g$, denoted
by $\mu(g)=E(Y_{i2}(g)|\mathcal{A}_{n})$, which captures the causal
relationship between peers’ and individual behavior as $g$ varies in
$\mathcal{G}$. We can also view $\mu(g)$ as an average dose-response function
that depends on the dose of treatment for the subset $\mathcal{A}_{n}$
(Forastiere et al., 2021). In the presence of network dependency, causal
estimands are commonly defined by conditioning on the specific network,
$\mathcal{A}_{n}$, without specifying the sampling mechanism from a larger
population. This is because describing the sampling mechanism is challenging
to conceptualize in network settings (Forastiere et al., 2022), and the
average of the potential outcome converges to a fixed quantity independent of
the particular network $\mathcal{A}_{n}$ under some conditions (Cristali and
Veitch, 2022). It is worth noting that the potential outcome can be defined
based on both $Y_{i1}$ and $G_{i}$, as in $Y_{i2}(Y_{i1}=y,G_{i}=g)$. However,
when the goal is to estimate peer influence while maintaining the observed
behaviors of node $i$ and varying only the behaviors of peers, the potential
outcome can be defined solely depending on $g$ (Egami and Tchetgen Tchetgen,
2023).

###  2.3 Identifying Assumptions

Given that the ACPI involves comparing two distinct quantities of potential
outcome at different levels of $g$, the identification of the average
potential outcome allows us to identify the ACPI. We outline the necessary
assumptions required for the identification of the average potential outcome
as following.

######  Assumption 1.

(Consistency)

Given a function of
$\phi:\mathbb{R}^{|\mathcal{N}_{i}|}\rightarrow\mathbb{R},\ \forall
i\in\mathcal{A}_{n}$,

$\displaystyle Y_{i2}=Y_{i2}\left(g\right)\text{ if
}g=\phi(\bm{Y}_{j1}:j\in\mathcal{N}_{i}).$ (1)

Assumption 1 implies that an individual’s potential outcome under the observed
peers’ behaviors at baseline is the outcome that will actually be observed for
that person.

######  Assumption 2.

(Neighborhood interference)

If
$\phi(\bm{Y}_{j1}:j\in\mathcal{N}_{i})=\phi(\bm{Y}^{\prime}_{j1}:j\in\mathcal{N}_{i})$
for all $i\in\mathcal{A}_{n}$, then

$\displaystyle
Y_{i2}\left(\phi(\bm{Y}_{j1})\right)=Y_{i2}\left(\phi(\bm{Y}^{\prime}_{j1})\right).$
(2)

Assumption 2 states that the potential outcomes take the same value when the
function values of different combinations of the peers’ behavior are equal.
Assumptions 1 and 2 can be interpreted as stable unit treatment on
neighborhood value assumption (SUTVA) in the presence of network interference
(Forastiere et al., 2021).

######  Assumption 3.

(Positivity)

$\displaystyle
0<P(G_{i}=g|Y_{i1}=y,\operatorname{\bm{X}}_{i}=x,U_{i}=u)<1\quad\text{for all
}g,u,x,\text{ and }y$ (3)

This assumption states that each individual has some chance of having a
possible value of $g\in\mathcal{G}$. However, this assumption of positivity
may not always be satisfied. For example, isolated nodes without any network
connection do not meet this assumption because the function values of $g$
cannot be calculated ($\mathcal{N}_{i}=\emptyset$). Isolated nodes could be
treated and analyzed separately or considered as $G_{i}=0$, as if they have
peers and their peers’ behaviors were 0 (Forastiere et al., 2021). This
scenario is beyond the scope of our paper.

######  Assumption 4.

(Latent unconfoundedness)

$\displaystyle
Y_{i2}(g){\perp\\!\\!\\!\perp}G_{i}|Y_{i1},U_{i},\operatorname{\bm{X}}_{i}\quad\forall
i,\ g\in\mathcal{G}.$ (4)

Assumption 4 implies that one’s potential outcome is independent of peers’
behaviors at baseline conditional on the individual behavior $Y_{i1}$,
observed confounder $\operatorname{\bm{X}}_{i}$ and unmeasured confounder
$U_{i}$. Therefore this assumptions suggests that conditioning on $U_{i}$ is
needed to control for confounding. However, as noted in Egami and Tchetgen
Tchetgen (2023), Assumption 4 while necessary, is not sufficient for potential
outcome identification because we do not observe $U_{i}$.

######  Assumption 5.

(Proxy variable for homophilous attributes)

There exists a proxy $Z_{i}$ for the unmeasured confounder $U_{i}$, and
$U_{i}$ is $Z_{i}$-measurable, satisfying the following condition:

$\displaystyle Y_{i2}\perp
A_{ij}|Z_{i},G_{i},Y_{i1},\operatorname{\bm{X}}_{i}.$ (5)

The Assumption 5 indicates that conditioning on $Z_{i}$ removes the effect of
conditioning on the $A_{ij}$, thereby preventing the opening of the backdoor
path. The measurability assumption implies that $U_{i}$ can be determined
using the information provided by $Z_{i}$ (Cristali and Veitch, 2022).

######  Theorem 1.

Under Assumptions 1-5, we have

$\displaystyle\operatorname{\mathbb{E}}\left[Y_{i2}(g)|\mathcal{A}_{n}\right]$
$\displaystyle=\operatorname{\mathbb{E}}\left[Y_{i2}|Y_{i1}=y,G_{i}=g,\operatorname{\bm{X}}_{i}=\operatorname{\bm{x}},Z_{i}=z\right]$
(6)
$\displaystyle=\operatorname{\mathbb{E}}\left[Y_{i2}|Y_{i1}=y,G_{i}=g,\operatorname{\bm{X}}_{i}=\operatorname{\bm{x}},Z_{i}^{*}=z^{*}\right]$
(7)

where $Z_{i}^{*}=WZ_{i}$ and $W$ is an orthogonal matrix satisfying
$W^{\prime}W=I$. Therefore, the average potential outcome can be identified by
leveraging a proxy variable $Z_{i}$ or a transformed proxy variable $Z^{*}$,
the observed data ($y$, $g$, $x$), and a set of nodes $i\in\mathcal{A}_{n}$.
The proof is provided in the Appendix.

The identifying assumptions we presented here, excluding the assumption on the
proxy variable $Z_{i}$ (Assumption 5), are similar to those in Egami and
Tchetgen Tchetgen (2023) and Liu and Tchetgen Tchetgen (2022), where an
observed variable associated with $U_{i}$ is selected for $Z_{i}$. However,
when $Z_{i}$ is unobserved, as it is in our case, we can regard a proxy
variable $Z_{i}$ as a latent variable that serves as a close approximation of
$U_{i}$. This proxy variable can be learned by applying a latent variable
social network model to estimate the latent structure that captures the
homophily (Louizos et al., 2017). Therefore, we further assume that we have
access to such model which allows us to sample from the posterior distribution
$f(Z_{i}|\mathcal{A}_{n})$ conditional on the observed network
$\mathcal{A}_{n}$.

There is an inherent trade-off between the bias that results from omitting the
unmeasured homophilous confounder $U_{i}$ and the bias from imprecisely
estimating or selecting the proxy variable $Z_{i}$. However, if the proxy
variable $Z_{i}$ contains sufficient information about $U_{i}$, accounting for
$Z_{i}$ can substantially mitigate bias compared to omitting it. Additionally,
as covariates measured in an observational study are not always eligible as
proxy variables for the unmeasured homophilous confounder, our emphasis in
this paper is on estimating $Z_{i}$ from the observed network. We leverage
existing latent variable network models (Dabbs et al., 2020) to estimate the
proxy variable that represents the underlying latent homophily in the network.
This estimated proxy variable enables the identification in Theorem 1
satisfying the assumptions presented in Section 2.3. In this paper, we focus
on a continuous latent variable $Z_{i}$, however the method can be extended to
accommodate discrete latent variables, such as block membership in a
stochastic block model (Fortunato, 2010; Holland et al., 1983; Karrer and
Newman, 2011).

##  3 Hierarchical Bayesian Model for the Estimation of Average Causal Peer
Influence Effects

In this section, we introduce estimation strategies for the causal estimand,
ACPI, within a Bayesian framework. To estimate the ACPI, we rely on the
posterior samples from $f(Z_{i}|\mathcal{A}_{n})$ utilizing a well established
latent variable network model, a latent space model with a binary tie value
(Hoff, 2005). We also demonstrate the eligibility of the resulting estimates
of latent space positions as proxy variables satisfying the required
conditions, and we extend our strategy to accommodate multiple network
settings. Further, after conditioning on the estimated $Z_{i}$, we incorporate
a Bayesian nonparametric model, for the outcome model. We select a Bayesian
additive regression tree (BART) for its ease-of-use and excellent predictive
performance (Chipman et al., 2010).

###  3.1 Latent Space Network Model

We consider the latent space model (LSM) (Hoff et al., 2002) which assumes
that the formation of network tie $A_{ij}$ from node $i$ to $j$ depends on
similarity between latent variables $Z_{i}$ and $Z_{j}$. In this setting, the
latent variable $Z_{i}$, is the position of node $i$ in a continuous metric
space (often the Euclidean metric), and the probability of forming a tie
$p(A_{ij}=1)$ is inversely related to the distance between $Z_{i}$ and
$Z_{j}$. Conditioning on these latent positions, which account for the
existing network structure such as homophily, transitivity and clustering,
network ties are assumed to be independent. The LSM can be written as

$\displaystyle A_{ij}\stackrel{{\scriptstyle
iid}}{{\sim}}\operatorname{Bernoulli}\left(p_{ij}\right),\quad\text{for
}i,j=1,\cdots,n$ (8)
$\displaystyle\operatorname{logit}\left(p_{ij}\right)=\alpha_{0}-\beta_{z}||Z_{i}-Z_{j}||_{2}$

where $Z_{i}\in\mathbb{R}^{d}$ is the position of individual $i$ in the latent
space. Denote $\bm{Z}=(Z_{1},\cdots,Z_{n})$. The intercept $\alpha_{0}$
represents a baseline probability of a tie for nodes in the same position and
$\beta_{z}$ is the corresponding slope coefficient for the distances between
$Z_{i}$ and $Z_{j}$. To estimate latent locations and parameters, both
likelihood-based inference and Bayesian inference are feasible approaches
(Hoff et al., 2002). In this paper, we will estimate LSM parameters within a
Bayesian framework, which allows for uncertainty quantification and is useful
for inference with small sample sizes (Lee and Song, 2004).

The prior distributions are generally specified as:

$\displaystyle Z_{i},Z_{j}\sim MVN_{d}(0,\sigma_{z}I)$ (9)
$\displaystyle\alpha_{0}\sim N(\mu_{\alpha},\sigma_{\alpha}^{2})$
$\displaystyle\beta_{z}\sim N(\mu_{\beta},\sigma_{\beta}^{2})$

where $Z_{i}$ and $Z_{j}$ are independently sampled from a multivariate normal
distribution with covariance matrix $\sigma_{z}I$ where $I$ is the
$d$-dimensional identity matrix. The intercept $\alpha_{0}$ and coefficient
$\beta_{z}$ are assigned normal prior distributions with hyperparameters
$(\mu_{\beta},\sigma_{\beta}^{2})$ and $(\mu_{\alpha},\sigma_{\alpha}^{2})$
for mean and variance, respectively.

The LSM can further be extended to model multiple networks using a
hierarchical latent space model (HLSM; Sweet et al. (2013)). By introducing a
network superscript $k$ to the variables in Equation (8), we define
$A_{ij}^{k}$ as the value of the tie from individual $i$ to individual $j$ in
network $k$. We denote the $k$-th adjacency matrix by $A^{k}$, where
$k=1,\cdots,K$, and $\bm{Z}^{k}=(Z_{1}^{k},\cdots,Z_{n_{k}}^{k})$ as the
latent positions for network $k$ with $n_{k}$ nodes. Following Sweet et al.
(2013), the HLSM can be written as

$\displaystyle A_{ij}^{k}\sim\operatorname{Bernoulli}\left(p_{ij}^{k}\right)$
(10)
$\displaystyle\operatorname{logit}\left(p_{ij}^{k}\right)=\alpha_{0}^{k}-\beta_{z}^{\prime}\left\|Z_{i}^{k}-Z_{j}^{k}\right\|.$

The ties are generated from the model based on each $\alpha_{0}^{k}$ and
$Z_{i}^{k}$, determining the structure of the kth network. To effectively
model networks that are independent yet similar, the following hierarchical
structure is considered:

$\displaystyle Z_{i}^{k},Z_{j}^{k}\stackrel{{\scriptstyle
iid}}{{\sim}}MVN_{d}(0,\sigma_{z}^{\prime}I)$ (11)
$\displaystyle\alpha_{0}^{k}\sim
N(\mu_{\alpha^{\prime}},\sigma_{\alpha^{\prime}}^{2})$
$\displaystyle\beta_{z}^{\prime}\sim
N(\mu_{\beta^{\prime}},\sigma_{\beta^{\prime}}^{2}).$

For $k=1,\ldots,K$, the parameters $\alpha_{0}^{k}$ and $\beta_{z}^{\prime}$
as well as network specific latent locations $Z_{i}^{k}$ are independently
drawn from a common population. This common population enables the effective
modeling of networks that are independent yet similar.

Both the LSM and the HLSM model tie formation as a function of the distance
between nodes. The distances between a collection of nodes in Euclidean space
are preserved under rotation, reflection, and translation (Hoff et al., 2002);
this implies that an infinite number of node positions results in the same
log-likelihood, given that
$\log\Pr(A|\bm{Z},\alpha_{0},\beta_{z})=\log\Pr(A|\bm{Z}^{*},\alpha_{0},\beta_{z})$
for any $\bm{Z}^{*}$ equivalent to $\bm{Z}$ through the transformation where
$\bm{Z}^{*}=(Z_{1}^{*},\cdots,Z_{n}^{*})$. Consequently, the latent positions
$\bm{Z}$ are identifiable only up to the class of distance preserving
transformations. However, this non-identifiability can be addressed by using a
Procrustes transformation (Borg and Groenen, 2005), to align each of the
latent position posterior draws with a fixed target position $\bm{Z}_{0}$ as a
post-processing step. The Procrustes transformation can be implemented by
targeting
$\bm{Z}^{*}=\arg\min_{W\bm{Z}}\operatorname{tr}\left(\bm{Z}_{0}-W\bm{Z}\right)^{\prime}\left(\bm{Z}_{0}-W\bm{Z}\right)$,
where $W$ is an orthogonal matrix $W\in\mathbb{R}^{d\times d}$ with
$W^{T}W=I$. This can easily be derived by
$\bm{Z}^{*}=\bm{Z}_{0}\bm{Z}^{\prime}(\bm{Z}\bm{Z}_{0}^{\prime}\bm{Z}_{0}\bm{Z}^{\prime})^{-1/2}\bm{Z}$
assuming $\bm{Z}_{0}$ and $\bm{Z}$ are both centered at the origin. The fixed
target positions $\bm{Z}_{0}$ can be set as the initial latent positions in
the MCMC algorithm by employing multidimensional scaling (MDS) on the observed
distance matrix. For the detail, we refer to Hoff et al. (2002); Adhikari et
al. (2021).

In an HLSM, the Procrustes transformation is performed on the posterior draws
of the latent positions $\bm{Z}^{K}$ for each network $k$ where each network
has a different fixed target positions $Z_{0}^{k}$ such that
$\bm{Z}^{*k}=\arg\min_{W^{k}\bm{Z}^{k}}\operatorname{tr}\left(\bm{Z}^{k}_{0}-W^{k}\bm{Z}^{k}\right)^{\prime}\left(\bm{Z}^{k}_{0}-W^{k}\bm{Z}^{k}\right)$
for $k=1,\cdots,K$. This is particularly important when $n_{k}\neq
n_{k^{\prime}}$ for $k,k^{\prime}\in(1,\cdots,K)$ and $k\neq k^{\prime}$. In
the multiple network setting, $\operatorname{\bm{U}}^{k}$ is identifiable as
$W^{k}Z^{k}$ and the rotation matrix $W^{k}$ varies across networks.

###  3.2 Outcome Model with Bayesian Additive Regression Trees

Let $\mathcal{D}_{i}=(Y_{i1},G_{i},Z_{i},\operatorname{\bm{X}}_{i})$ denote a
vector of predictors for node $i$, where $Z_{i}$ represents the latent
location sampled from the LSM. Then the outcome model for behavior at follow-
up, $Y_{i2}$ for node $i\in(1,\cdots,n)$ is written as

$Y_{i2}=f(\mathcal{D}_{i})+\epsilon_{i}$

where $\epsilon_{i}\sim\mathrm{N}(0,\tau^{2})$. To model the unknown mean
function $f(\cdot)$, we consider BART, which offers a flexible and
nonparametric approach to model multidimensional regression bases being
learned directly from the observed data(Chipman et al., 2010). Such a
nonparameteric approach accommdates a variety of nonlinear associations
between peer influence and behavior.

The BART framework introduced by models $f(\mathcal{D}_{i})$ as a sum of $m$
trees

$\displaystyle
f(\mathcal{D}_{i})=\sum_{t=1}^{m}g(\mathcal{D}_{i};\mathcal{T}_{t},\mathcal{M}_{t})\
,$ (12)

where $\mathcal{T}_{t}$ is a binary tree structure with $n_{t}$ leaf nodes and
$\mathcal{M}_{t}=\\{\mu_{t1},\ldots,\mu_{tn_{t}}\\}$ is a collection of
parameters for $n_{t}$ leaves of the $t$th tree $\mathcal{T}_{t}$. The
function $g(\mathcal{D}_{i};\mathcal{T}_{t},\mathcal{M}_{t})$ returns
$\sum_{\ell=1}^{n_{t}}\mu_{t\ell}\ \eta(\mathcal{D}_{i};\mathcal{T}_{t},\ell)$
where $\eta(\mathcal{D}_{i};\mathcal{T}_{t},\ell)$ is the indicator that
$\mathcal{D}_{i}$ is associated with leaf node $\ell$ in tree
$\mathcal{T}_{t}$, that is,
$g(\mathcal{D}_{i};\mathcal{T}_{t},\mathcal{M}_{t})=\mu_{t\ell}$. The tree
$\mathcal{T}_{t}$ is constructed using interior splitting rules of the form
$[d_{j}\leq C_{b}]$ at each branch node $b$, inducing a partition on the
covariate space. The cutpoint $C_{b}$ of the splitting rule is assigned a
uniform prior over its possible values. To avoid overfitting, each tree is
given a _regularization prior_ which shrinks the individual leaf parameters to
$0$, ensuring that each tree contributes only a small portion of the overall
fit. The regularization is governed by the depth of each tree and the prior
distribution of leaf parameters. The probability that each node at depth $d$
is non-terminal is given by $\gamma(1+d)^{-\beta}$ for hyperparameters
$\gamma\in(0,1)$ and $\beta\in[0,\infty)$. The default prior specification
(Chipman et al., 2010) sets $\gamma=0.95$ and $\beta=2$, which encourages
shallow trees, rarely exceeding depth $2$–$3$. The leaf parameters
$\mu_{t\ell}$ are assigned iid $N(0,\sigma^{2}_{\mu}/m)$ priors, ensuring that
the prior variance of $f(\bm{x})$ remains constant as a function of $m$.

In extending to the multiple network setting, while considering varying
$W^{k}$ across networks in the outcome model, we consider a semi-parametric
BART model given as

$\displaystyle Y^{k}_{i2}$
$\displaystyle=\sum_{t=1}^{m}g(\mathcal{D}^{k}_{i};\mathcal{T}_{t},\mathcal{M}_{t})+(\beta_{z}^{k})^{T}\bm{Z}_{i}^{k}+\epsilon_{i}^{k}$
(13)

where
$\mathcal{D}^{k}_{i}=(Y_{i1}^{k},G_{i}^{k},\operatorname{\bm{X}}_{i}^{k})$,
$\beta_{z}^{k}\in\mathbb{R}^{d}$ and $\epsilon_{i}^{k}\sim N(0,\tau^{\prime
2})$. The coefficient $(\beta_{z}^{k})$ varies with network $k$ and adapts to
the rotation matrix $W^{k}$ such that
$(\beta_{z}^{k})^{T}\bm{Z}_{i}^{k}=(\beta_{z})^{T}W^{k}\bm{Z}_{i}^{k}$. By
using a semi-parametric BART model, we account for network-specific
homophilous attributes which are not shared across the networks via the random
effect $(\beta_{z}^{k})$. This model adopts a similar strategy to mixedBART
(Spanbauer and Sparapani, 2021), where a parametric random effect term
$(\beta_{z}^{k})^{T}\bm{Z}_{i}^{k}$ is introduced into the BART model
analogous to random effect term in a linear mixed effect model. The
$\sum_{t=1}^{m}g(\mathcal{D}^{k}_{i};\mathcal{T}_{t},\mathcal{M}_{t})$
represents the population fixed effect regression relationship, shared by all
of the observations, analogous to a fixed effect term in a linear mixed effect
model. Therefore, this semi-parametric BART model effectively captures the
nonlinear relationship between peer influence and individual behavior, which
is shared by all observations, while accounting for varying rotation matrix
$W^{k}$ across networks. For the parameter $\beta_{z}^{k}$, a multivariate
normal prior distribution $\beta_{z}\sim N_{d}(0,\sigma_{\beta}I)$ is used.
The details of our Bayesian computation and estimation algorithm are provided
in Appendix C.

Algorithm 1 LSM-BART

1:for $j=1$ to $J$ do

2: (Latent location stage)

3: Draw a sample of $\alpha_{0}^{j}\sim
p(\cdot|\mathbf{Z}^{j-1},\beta_{z}^{j-1},A)$ where
$\mathbf{Z}^{j-1}=(Z_{1}^{j-1},\cdots,Z_{n}^{j-1})$

4: Draw a sample of $\beta_{z}^{j}\sim
p(\cdot|\mathbf{Z}^{j-1},\alpha_{0}^{j},A)$.

5: for $i=1$ to $N$ do

6: Update $Z_{i}^{j-1}$ given $Z_{i}^{j},\alpha_{0}^{j},\beta_{z}^{j}$ and $A$

7: end for

8: (Outcome stage)

9: for $m=1$ to $M$ do

10: Set $R^{j}_{mi}=y_{i2}-\sum_{q\neq
m}g\left(D_{i};\mathcal{T}_{q},\mathcal{M}_{q}\right)$ for $i=1,\ldots,n$
where $D_{i}=(Y_{i1},G_{i},\mathbf{Z}_{i}^{j}$).

11: Propose a new tree $T_{m}^{*}\sim
Q\left(\mathcal{T}_{m}\rightarrow\mathcal{T}_{m}^{*}\right)$.

12: Compute the acceptance ratio  $\displaystyle
r=\frac{\mathscr{L}\left(\mathcal{T}_{m}^{*}\right)Q\left(\mathcal{T}_{m}^{*}\rightarrow\mathcal{T}_{m}\right)}{\mathscr{L}\left(\mathcal{T}_{m}\right)Q\left(\mathcal{T}\rightarrow\mathcal{T}_{m}^{*}\right)},$
where
$\mathscr{L}\left(\mathcal{T}_{m}\right)=p\left(\mathcal{T}_{m}\mid\bm{R}_{m},\sigma\right)\times
p\left(\mathcal{T}_{m}\right)$ and
$\bm{R}_{m}=(R^{j}_{m1},\ldots,R^{j}_{mn})$.

13: Set $\mathcal{T}_{m}=\mathcal{T}_{m}^{*}$ if $U\leq\min(1,r)$ where
$U\sim\operatorname{Uniform}(0,1)$, otherwise retain $\mathcal{T}_{m}$.

14: Sample the $\mu_{t\ell}$’s from their full conditional distributions.

15: end for

16: Update $\tau$ from its full conditional distribution.

17: (Monte Carlo estimation stage)

18: for $g\in\mathcal{G}$ do

19: Predict $\hat{Y}_{i2}(g)$ given $Z_{i}$,
$\bm{\mathcal{T}}=(\mathcal{T}_{1},\ldots,\mathcal{T}_{m}),$ and
$\bm{\mathcal{M}}=(\mathcal{M}_{1},\ldots,\mathcal{M}_{m})$

20: end for

21: Compute average the potential outcomes over all units as
$\hat{\mu}(g)=\frac{1}{N}\sum_{i=1}^{N}\hat{Y}_{i2}(g)$ for all
$g\in\mathcal{G}$

22: Compute ACPI as $\tau(g,g^{\prime})=\hat{\mu}(g)-\hat{\mu}(g^{\prime})$
for all $g,g^{\prime}\in\mathcal{G}$

23:end for

###  3.3 Two-step Estimation

We estimate our model using a Markov chain Monte Carlo (MCMC) algorithm
consisting of two steps: 1) the update of latent locations and related LSM
parameters, 2) the update of the parameters in the outcome model. Let
$\bm{\theta}_{z}$ and $\bm{\theta}_{y}$ be the vectors of parameters of
network (LSM) and the outcome model, respectively. The full conditional
distributions for each parameter are

$\displaystyle
p\left(\bm{\theta}_{y}\mid\bm{Y}_{2},\bm{Y}_{1},\mathbf{X},\bm{\theta}_{z},A\right)\propto
p\left(\bm{Y}_{1},\bm{Y}_{2},\mathbf{Z},\mathbf{X},A\mid\mathbf{Z},\bm{\theta}_{z},\bm{\theta}_{y}\right)p\left(\bm{\theta}_{y}\right)$
$\displaystyle
p\left(\mathbf{Z},\bm{\theta}_{z}\mid\bm{Y}_{2},\bm{Y}_{1},\mathbf{X},\bm{\theta}_{y},A\right)\propto
p\left(\bm{Y}_{1},\bm{Y}_{2},\mathbf{Z},\mathbf{X},A\mid\mathbf{Z},\bm{\theta}_{z},\bm{\theta}_{y}\right)p\left(\mathbf{Z},\bm{\theta}_{z}\right)$
(14)

respectively. Within the Bayesian framework, the joint posterior distribution
of $\mathbf{Z}$ and $\bm{\theta}_{z}$ is in part informed by the update in
outcome model. This may raise concerns when 1) connections between individuals
were established before the behaviors are disseminated, and 2) model
misspecification in the relationship between outcome and latent locations
results in incorrect causal estimates. This is similar to model feedback
(Hoshino, 2008; Zigler et al., 2013) where the posterior distribution of the
parameters in the propensity score model is partially informed by the outcome
model.

To address these concerns, we replace the full conditional distribution
$\eqref{eq:posterior}$ with an approximate conditional distribution given as

$\displaystyle
p\left(\mathbf{Z},\bm{\theta}_{z}\mid\bm{Y}_{2},\bm{Y}_{1},\mathbf{X},\bm{\theta}_{y},A\right)\propto
p\left(A\mid\mathbf{Z},\bm{\theta}_{z}\right)p\left(\mathbf{Z},\bm{\theta}_{z}\right)$

employing a similar strategy as in Hoshino (2008); Kaplan and Chen (2012);
McCandless et al. (2010); Bayarri et al. (2009).

Samples are drawn from the posterior predictive distribution of the potential
outcomes $Y_{i2}(g)$ for each $i$. Subsequently, for each sample, we compute
the average potential outcome $\mu(g)$ by averaging $Y_{i2}(g)$ across all
nodes. The ACPI is then computed by comparing the average potential outcomes
at different levels. A description of the estimation process using the LSM is
provided in Algorithm 1, and the algorithm with the HLSM is available in the
Appendix.

##  4 Simulation Study

To assess the performance of our proposed framework for estimating the ACPI,
we conduct a series of three simulation studies where we compared our method
to a naive estimator that does not adjust for latent homophily and an oracle
estimator assuming that the latent homophilous confounder $U_{i}$ is observed.
In the first simulation setting in Section 4.1, we assess how the strength of
network homophilly, as measured by the association between the latent
attribute and the likelihood of tie formation, impacts ACPI estimation. In the
second set of simulations in Section 4.2, we investigate how the strength of
unmeasured confounding, quantified by the dependence between the latent
attribute and the baseline behavior, may impact the estimation. The first two
simulation settings focus on single network data. In the last simulation
setting in Section 4.3, we investigate ACPI estimation with multiple networks,
similar to our motivating example.

We first present our data generating process for a single network using the
LSM. To incorporate homophily through clusters, we generate attributes
$\operatorname{\bm{U}}_{i}\in\mathbb{R}^{2}$ for each node $i$ using Gaussian
mixture model with $J$ components, such that the $j$th component has a mean of
$\mu_{j}$ and a diagonal variance-covariance matrix of $\sigma^{2}_{j}I$, as
$\bm{U}_{i}\sim\sum_{j=1}^{J}\phi_{j}MVN_{2}(\mu_{j},\sigma^{2}_{j}I),$ where
$\phi_{j}$ is the mixture component weight subject to the constraint
$\sum_{j=1}^{J}\phi_{j}=1$. In all simulations, we allocate each node to a
component by sampling the corresponding component $j$ using a
$Multinomial(n;p_{1},\cdots,p_{J})$, where $p_{j}=1/J$. We fix $J=4$ and the
mean vector for each component of the Gaussian mixture model as
$\mu_{1}=(-1.5,-1.5)$, $\mu_{2}=(1.5,-1.5)$, $\mu_{3}=(-1.5,1.5)$, and
$\mu_{4}=(1.5,1.5)$. We set $\phi_{j}$ to $1/J$ to assign equal weight to each
component. Next we generate an undirected network following LSM as specified
in (8) assuming $U_{i}=Z_{i}$ with $\alpha_{0}=6$ and $\beta_{z}=5$. In all
simulations, we fix $\beta_{z}$ as we demonstrate variations in homophily
strength by varying $\sigma_{j}$ instead. Further, for all simulations with a
single network we fix the number of nodes $n$ to 100. To generate multiple
networks for simulations in Section 4.3, we use the HLSM in (10) by setting
$\alpha_{0}^{k}=3$ in (10). Each network is generated following the same
strategy as in the single network setting using a Gaussian mixture model with
$n_{k}=40$ and $J=4$ where the number of nodes within each network was chosen
to mimic our motivating example.

Given our 100-node networks, we generate individual behaviors at baseline as
$Y_{i1}=5+\beta_{u_{1}}U_{i1}+\beta_{u_{2}}U_{i2}+\epsilon_{i1}$, where
$\epsilon_{i}\sim N(0,1)$. We further specify a nonlinear outcome model to
generate individual behaviors at follow-up as
$Y_{i2}=10+3Y_{i1}+10\text{sigmoid}(10(G_{i}-5))+\beta^{\prime}_{u_{1}}U_{i1}+\beta^{\prime}_{u_{2}}U_{i2}+\epsilon_{i2}$,
where $\text{sigmoid}(x)=1/(1+e^{-x})$. With these true underlying values, the
population average potential outcome is
$\operatorname{\mathbb{E}}(Y_{i2}(g))=25+10\text{sigmoid}(10(g-5))$. Details
of the data generating process in the multiple networks setting are provided
in Appendix D

To investigate the influence of homophily on ACPI, we vary $\sigma_{j}$ in the
Gaussian mixture model for $U_{i}$ in Section 4.1, whereas we fix
$\sigma_{j}=0.3$ in the other two simulation settings. Similarly, to assess
the impact of varying strength of association of $U_{i}$ with $Y_{i1}$ and
$Y_{i2}$ (simulations in Section 4.2), we vary $\beta_{u_{1}}$,
$\beta_{u_{2}}$, $\beta^{\prime}_{u_{1}}$ and $\beta^{\prime}_{u_{2}}$. For
other simulations they are fixed at
$\beta_{\mathbf{u}}=(\beta_{u_{1}},\beta_{u_{2}})=(0.5,0.5)$ and
$\beta^{\prime}_{\mathbf{u}}=(\beta^{\prime}_{u_{1}},\beta^{\prime}_{u_{2}})=(10,10)$.
Finally, we vary the number of networks in Section 4.3.

To sample from the posterior distribution of the latent locations $Z_{i}$
using the LSM, we set the hyperparameters as $\sigma_{z}=5,\ \mu_{\alpha}=0,\
\mu_{\beta}=0,\ \sigma_{\alpha}=10,\ \sigma_{\beta}=10$. For the multiple
networks setting, we set $\sigma_{z}^{\prime}=5,\ \mu_{\alpha^{\prime}}=0,\
\mu_{\beta^{\prime}}=0,\ \sigma_{\alpha^{\prime}}=10,\
\sigma_{\beta^{\prime}}=10$. We next estimate the posterior mean (with 95%
posterior intervals) of the ACPI using either BART (for single network) or
Semi-BART (for multiple networks) with 10000 MCMC iterations with 5000 burn-
in. We set the number of trees $m$ to 200. While it is possible to include the
additional confounders $X_{i}$ that are independent of network tie formation
yet associated with $Y_{i1}$ and $Y_{i2}$, our simulation excludes $X_{i}$ to
precisely investigate how unobserved homophily can bias ACPI estimation and
whether these latent homophilous attributes can be estimated by conditioning
on the observed network data. Further details on the estimation can be found
in Appendix D. We use the _HLSM_ R package to estimate latent locations.

The average potential outcome $\mu(g)$ is calculated by averaging $Y_{i2}(g)$
across all nodes, as outlined in Algorithm 1. This computed $\mu(g)$, is also
examined in the simulation studies. It is important to note that while
$\mu(g)$ can be computed for all $g\in\mathcal{G}$, the primary estimand of
interest, ACPI, involves comparisons between the average potential outcomes
$\mu(g)$ of two distinct treatment levels, $g$ and $g^{\prime}$. To evaluate
the ACPI estimation with two levels, we partition the range of interest into
10 equally spaced intervals, represented as $(g_{1},\ldots,g_{10})$. We then
compute the ACPI, defined as $\tau(g_{q+1},g_{q})=\mu(g_{q+1})-\mu(g_{q})$,
for each consecutive pair $q=1,\ldots,9$. We set this range for peer behaviors
from 4 to 6, such as $(g_{1},\ldots,g_{10})=(4,4.22,\ldots,5.77,6)$, which
covers the 95% range of observed peer behaviors in all simulation settings.

For each data generating model, we replicate the simulation 10 times, and
three different criteria for evaluation are computed at each replication. The
first criteria is the root mean squared error (RMSE) in estimating the ACPI
with the posterior mean. The second criteria is the coverage or the proportion
of 95% credible intervals that capture the true value of the ACPI. The third
criteria is the average length of the 95% credible interval for the ACPI.
Then, our evaluation criteria are averaged over 10 replications and the entire
grid. Specifically,
$RMSE=\sqrt{\sum_{r}\sum_{q}\left(\hat{\tau}^{r}(g_{q+1},g_{q})-\tau^{r}(g_{q+1},g_{q})\right)^{2}/RQ}$,
where $R$ is the number of replications and $Q$ is the number of grid points
in the range $(Q=10)$. The $\hat{\tau}^{r}(\cdot)$ represents the posterior
mean at the $r$th replicate. The coverage and length are defined in a manner
similar to RMSE.

Figure 3: Simulation 1: The generated networks from data generating process
with $\sigma_{j}\in(0.1,0.4,0.7,1)$.

###  4.1 The Strength of Homophily and ACPI

Figure 4: Simulation 1: The RMSE, coverage, and length of ACPI estimation from
oracle, LSM, and naive estimator are evaluated based on 10 replications of the
simulation at each standard deviation value.

In this simulation study, we examine the impact of the homophily’s strength on
ACPI estimation, with the hypothesis that strong clustering should induce
strong homophily as well as increased oppotunities for influence (Warren et
al., 2020; Malik et al., 2016). When we generate latent locations, the
variance of latent locations within the same component is controlled by
$\sigma_{j}$, which can be interpreted as the strength of the homophily. We
evaluate the impact of these different homophily strengths within a component
by varying $\sigma_{j}\in(0.1,0.2,\ldots,0.9,1)$. Examples of generated
networks with $\sigma_{j}=(0.1,0.4,0.7,1)$ are illustrated in Figure 5. When
$\sigma_{j}=0.1$, nodes within the same component are more likely to be
connected while connections between components are limited. Additionally,
smaller variations or closer proximity between nodes within a component imply
that the nodes share more similar attributes. In contrast, at $\sigma_{j}=1$,
nodes still tend to connect within their component, maintaining a homophilous
network, but the distinction between components becomes less apparent, leading
to increased connections between components.

The average RMSE, coverage, and interval length of the ACPI estimation
performance is represented in Figure 4. Both the LSM and Oracle exhibit
similar performances in terms of RMSE and coverage, with the LSM achieving
significantly narrower interval lengths compared to the Naive estimator. This
suggests that, regardless of the strength of homophily, including latent
locations leads to an improvement in mitigating bias in the estimation of ACPI
while also narrowing the ACPI credible intervals. Another observation is that
both Naive and LSM estimators perform poorly when standard deviations are
small ($\sigma_{j}\leq 0.3$), and the performance is stabilized for
($\sigma_{j}>0.3$) across all criteria.

Figure 5: Simulation 1: The scatter plot of observed individual behavior at
follow-up ($y_{i2}$) versus the (observed) average of peers’ behaviors with
$\sigma_{j}\in(0.1,0.4,0.7,1)$. Figure 6: Simulation 1: The estimated average
potential outcome $\mu(g)$ from oracle, LSM, and naive estimator with standard
deviations $\sigma_{j}$ of 0.1, 0.4, 0.7, and 1.0. The 95% confidence interval
is also represented by shading along with the true average potential outcome
(dotted line) obtained from 10 replications.

The underperformance with smaller standard deviations is primarily attributed
not to estimation inaccuracies but to a violation of the positivity
assumption. To see this violation, Figure 5 presents the observed individual
behaviors colored by each component, and Figure 6 shows the estimated average
potential outcomes with $\sigma_{j}$ = (0.1, 0.4, 0.7, 1). Given
$\sigma_{j}=0.1$, the observed behaviors in component 1 are less likely to be
similar to those in component 2 in Figure 5. That is, in the presence of
strong homophily, behaviors are more distinct across different components in
Gaussian mixture model as baseline behaviors are determined by the homophilous
attributes $\mathbf{U}_{i}$. Also, there exist ranges in which peers’
behaviors are not observed leading to biases in average potential outcome
estimation. It is important to note that behaviors observed from each
component do not cover the entire range of $g$. However, in Figure 6, our
oracle and LSM estimators effectively recover true average potential outcome,
except for strong homophily such as $\sigma_{j}=0.1$. Conversely, the
estimated average potential outcome from Naive estimator remains constant
irrespective of the strength of homophily, exhibiting the largest bias. This
consistent performance arises due to the Naive estimator’s ignorance of the
homophily effect.

###  4.2 Coefficient of Unmeasured Confounding and ACPI

Figure 7: Simulation 2: The estimated average potential outcome from oracle,
LSM, and naive estimator varying
$\beta_{\mathbf{u}}^{\prime}\in(1,\ldots,10)$, represented by different
colors. The population average potential outcome $\mu(g)$ is outlined with a
dotted line. Results for $\beta_{\mathbf{u}}=0$ are shown in the first row,
and for $\beta_{\mathbf{u}}=0.5$ in the second row. The 95% confidence
interval is also represented by shading along with the true average potential
outcome (dotted line) obtained from 10 replications. Figure 8: Simulation 2:
The RMSE, coverage, and length of ACPI estimation from oracle, LSM, and naive
estimator are evaluated based on 10 replications of the simulation at each
coefficient value $\beta^{\prime}_{\mathbf{u}}$. Results are shown with
$\beta_{\mathbf{u}}=0$ (top) and $\beta_{\mathbf{u}}=0.5$ (bottom).

In our second simulation study, we explore impact of the effect size of latent
homophilous attributes in behaviors at baseline and follow-up on the
estimation of ACPI by varying coefficients of $U_{i}$ in outcome models for
$Y_{i1}$ and $Y_{i2}$. We consider two scenarios: i) $U_{i}$ has no influence
on $Y_{i1}$ with $\beta_{\mathbf{u}}=(\beta_{u_{1}},\beta_{u_{2}})=(0,0)$ and
$\epsilon_{i2}\sim N(0,2^{2})$, and ii) $U_{i}$ has influence on $Y_{i1}$ with
$\beta_{\mathbf{u}}=(\beta_{u_{1}},\beta_{u_{2}})=(0.5,0.5)$ and
$\epsilon_{i1}\sim N(0,1)$. The first scenario indicates that the latent
homophilous attributes $\mathbf{U}_{i}$ are independent of $Y_{i1}$ implying
no arrow from from $U_{i}$ to $Y_{i1}$ in Figure 2 while keeping an arrow from
$U_{i}$ to $Y_{i2}$. This implies that behaviors at baseline are randomly
assigned, but the network ties still depend on latent homophilous attributes.
The standard deviation of $\epsilon_{i2}$ in the first scenario is increased
to ensure that the range of $Y_{i1}$ and subsequently $G_{i}$ is similar to
their range in the second scenario. For both scenarios, we also vary the value
of $\beta^{\prime}_{u_{1}}=\beta^{\prime}_{u_{2}}\in(1,2,\cdots,9,10)$, to
investigate the influence of the effect size of the latent homophilous
attribute in the behaviors at baseline on the estimation of average potential
outcome as well as ACPI.

The results of average potential outcome estimation are presented in Figure 7
for the first scenario (top) and the second scenario (bottom). In the first
scenario where $Y_{i1}$ are randomly assigned, the magnitude of the
coefficient $\beta^{\prime}_{u}$ in the follow-up outcome model does not have
any effect on the average potential outcome estimation across all three
methods. Also, the estimated average potential outcomes from all methods
capture ground truth well. However, the confidence interval obtained from 10
replications in the Naive estimator, is significantly wider, especially as
$\beta^{\prime}_{u}$ increases. This observation suggests that ignoring the
latent homophilous attribute results in a wider confidence interval when
$\operatorname{\bm{U}}_{i}{\perp\\!\\!\\!\perp}Y_{i1}$, thereby indicating
increased uncertainty in the estimation. In this result, ignoring latent
attributes does not contribute to significant bias in the estimation of the
average potential outcome when peers’ behaviors at follow-up are assigned
randomly; instead, it merely increases the variance. On the other hand, in the
second scenario (bottom panel in Figure 7), where the latent homophilous
attributes $\operatorname{\bm{U}}_{i}$ is involved in the behavior at
baseline, Naive estimator shows that the the bias in estimating average
potential outcome increase as value of the $\beta^{\prime}_{\mathbf{u}}$
increases. However, our proposed LSM estimator is effective in reducing this
type of bias by adjusting for the latent locations. Although the bias in the
estimated average potential outcome slightly increases as
$\beta^{\prime}_{\mathbf{u}}$ increases, it remains markedly lower compared to
the Naive estimator. Moreover, the performance of LSM estimator is comparable
to the result of Oracle estimator. This result demonstrates that adjusting for
the latent locations $Z_{i}$ effectively mitigates bias caused by homophily,
especially when the baseline behavior is also impacted by latent homophilous
attributes.

Next, we compute the average RMSE, coverage, and interval length, represented
for ACPI estimation in Figure 8 for the first scenario (top) and second
scenario (bottom). For both scenarios, the RMSE from the Naive estimator
strictly increases as $\beta^{\prime}_{\mathbf{u}}$ increases, while the
average coverage is consistently lower than that of its competitors.
Additionally, the CI length from Naive estimator is consistently larger than
those of both the LSM and the Oracle method,regardless of the values of
$\beta^{\prime}_{\mathbf{u}}$. Especially, RMSE from LSM estimator is
significantly reduced compared to that from Naive estimator.

###  4.3 Number of Networks and ACPI: Extension to HLSM

In this section, we extend our approach to the multiple networks setting and
investigate the impact of the number of networks on the estimation of ACPI. We
varied the number of networks $k\in(2,4,\ldots,12,14)$ and generated network
data from the HLSM given in (10).

Figure 9: Simulation 3: The RMSE, coverage and length of ACPI are evaluated
based on 10 replications of the simulation varying the number of networks with
$k\in(2,4,\ldots,12,14)$.

We compute the average RMSE, coverage, and interval length for each $k$,
represented in Figure 9. For all evaluated criteria, the Naive estimator
demonstrates inferior performance, exhibiting the highest RMSE and lowest
coverage. In contrast, the LSM-semiBART estimator significantly reduces the
RMSE and provides improved coverage. This indicates that incorporating
adjustments for latent locations effectively reduces bias associated with
homophily and decreases uncertainty within a multi-network context. Another
finding is that for all estimators, both RMSE and length remain relatively
consistent when $K\geq 6$. This indicates that a sample comprising 6 networks
($N=n_{k}\times K=40\times 6=240$) is adequately large to ensure robust
estimation of the ACPI in our simulation setting.

Unlike the simulations in single network settings described in Sections 4.1
and 4.2, the coverage of the LSM-semiBART estimator is smaller compared to
that of the Oracle estimator. This decrease may be attributed to the expansion
of the parameter space. For a fixed set of nodes, as the number of networks
($k$) increases, the dimension of the parameter
$\beta^{k}_{\mathbf{z}}=(\beta^{k}_{z_{1}},\beta^{k}_{z_{2}})$ also increases,
accounting for variations in transformations between networks. This expanded
parameter space likely results in a decrease in coverage as $k$ increases.
Nevertheless, the LSM-semiBART estimator significantly outperforms the Naive
estimator across all evaluation criteria, and the RMSE remains overall
consistent, regardless of the number of networks $k$. This performance
indicates that our proposed LSM-semiBART model effectively estimates the
shared ACPI across networks by accounting for the transformed latent locations
across different networks.

##  5 Real Data Analysis

To investigate the causal peer influence effect on beliefs about mathematics
using advice-seeking networks, we apply our proposed method to estimate the
ACPI. Recall from Section 1.2, our data include 14 advice-seeking networks
among elementary school staff in a school district. We hypothesize that
teachers’ beliefs in 2013 were shaped not only by their own beliefs but also
by those of their peers in 2012. For this analysis, from a total of 14
schools, we selected 11 schools with more than 20 respondents (nodes) and
assumed an undirected relationship, where all ties are considered to be
reciprocated. We use the advice-seeking network data collected in Spring 2013,
as teachers were asked to nominate others based on their interactions over the
past year.

Following the LSM-semiBART approach for multiple networks to assess the
average potential outcome $\mu(g)$, we define the range for $g$ to be from
$15.6$ to $22.6$, which empirically encompasses the 95% range of observed
average peer belief scores. Then we partition this range into 10 equally
spaced intervals such as
$(g_{1},g_{2},\cdots,g_{9},g_{10})=(15.6,16.0,\cdots,22.2,22.6)$ and evaluate
the average potential outcome $\mu(g_{q})$ for $q=1,\ldots,10$ and ACPI
$\tau(g_{q+1},g_{q})=\mu(g_{q+1})-\mu(g_{q})$ for $q=1,\ldots,9$.

The left panel of Figure 10 illustrates the estimated average belief scores
(simply referred to as average belief scores) in 2013 versus the average
peers’ belief score $g$ in 2012, estimated from the Naive estimator and LSM-
semiBART. These average belief scores are calculated across all schools, with
detailed estimates for each school provided in Figure 14 in the appendix. The
Naive estimator exhibits a fluctuating pattern in the average belief scores,
along with a wider 95% credible interval, compared to the LSM-semiBART
estimator. This pattern implies that neglecting the homophily effect
introduces bias and increases uncertainty in the estimating average belief
scores.

The right panel of Figure 10 represents the estimated ACPI
$\tau(g_{q+1},g_{q})$ for $q=1,\cdots,9$ from Naive and LSM-semiBART
estimators. The Naive estimator leads to ACPI being either overestimated or
underestimated, with fluctuations. Additionally, neglecting the latent
homophily results in wider credible intervals, thus diminishing the
reliability of the estimation. Although the impact of latent homophily on
estimating the ACPI in our dataset is not markedly significant, it helps
mitigate the bias introduced by overlooking the latent homophily and enhances
the precision of our ACPI estimates.

Figure 10: The estimated average belief scores in 2013 across all school from
the Naive and LSM-semiBART estimator is represented with 95% credible interval
(left). The estimated ACPI is depicted with 95% credible interval from each
estimator (right). The rug plot shows the observed distribution of peers’
belief score.

Finally to ensure that the HLSM specification provides a reasonable
approximation of the latent homophilous attribute, we perform goodness of fit
checks using the posterior predictive distributions of the networks. The
visualization of the network summaries (tie density and transitivity) of the
predicted networks in Figure 16 demonstrates that the network ties from the
prediced networks closely resemble those of the observed networks.

In our analysis of real data, we observed a lack of significant ACPI. Previous
studies, which employed different models and utilized latent locations in
various ways with directed networks Spillane et al. (2018); Sweet and Adhikari
(2020), also reported minimal evidence of a peer influence effect on student-
centered beliefs on average across schools. This insignificance of ACPI may be
attributed to the limited duration of the observation period.

##  6 Discussion

Social network structure is caused by either observed or unobserved homophily,
which can result in biased estimates of causal peer influence. Our proposed
LSM-BART (LSM-semiBART) estimator demonstrated that this bias can be mitigated
by adjusting for the latent locations of each node. We employed BART to
estimate the nonlinear relationship between individual behavior and peers’
behaviors and showed how our framework can be extended from a single network
to a multiple networks setting. Our estimation process leverages a two-step
Bayesian procedure, which enables the quantification of uncertainty in both
latent location and potential outcome estimation. Through a series of
simulation studies, we found that neglecting latent homophily led to biased
estimates of the ACPI, whereas adjusting for latent locations effectively
mitigated this bias. In the context of multiple networks, the semi-parametric
BART model effectively explained the variation in transformations across
networks as well as the nonlinear relationship between individual behavior and
peers’ behaviors.

In our real data analysis, we estimated the ACPI to gain insights into changes
in teachers’ beliefs about mathematics that resulted from peer influence.
Although the impact of latent homophily on the ACPI was not significant, we
observed that overlooking latent homophily can lead to either underestimation
or overestimation of the ACPI, accompanied by considerable uncertainty. Given
the potential benefits of accurately identifying the ACPI for enhancing
teaching quality, these findings highlight the importance of incorporating
latent locations to better understand homophily among teachers and applying
appropriate regression models. While node-level covariates could potentially
account for part of the latent homophily, we focus more on estimating latent
locations using solely the observed network to account for latent homophily.

Our framework focuses on identifying the effects of peer behaviors on
individual behaviors, where both are measured on a continuous scale. This
approach represents a generalized version of the study on causal spillover
effects in social networks (Forastiere et al., 2021, 2022) by transitioning
from a binary to a continuous treatment covariate. Thus, our framework can be
applied to spillover effect studies, especially when latent homophily exists
and the relationship between the outcome and spillover effect is nonlinear.
For future work, our framework can be extended to identify heterogeneous
causal peer influence effects using BART. While our current focus is on the
nonlinearity between individual behaviors and their peers’ behaviors, the
complex interactions among peer influence, confounders, covariates, and
treatment can be identified by applying BART.

##  7 Acknowledgement

We thank Jim Spillane and the Distributed Leadership Studies at Northwestern
University http://www.distributedleadership.org, including the NebraskaMATH
study, for the data used in this paper.

## References

  * Napier et al. [1984] Ted L Napier, Richard Goe, and Douglas C Bachtel.  An assessment of the influence of peer association and identification on drug use among rural high school students.  _Journal of Drug Education_ , 14(3):227–248, 1984\.
  * Henneberger et al. [2021] Angela K Henneberger, Dawnsha R Mushonga, and Alison M Preston.  Peer influence and adolescent substance use: A systematic review of dynamic social network research.  _Adolescent Research Review_ , 6(1):57–73, 2021\.
  * Christakis and Fowler [2007] Nicholas A Christakis and James H Fowler.  The spread of obesity in a large social network over 32 years.  _New England journal of medicine_ , 357(4):370–379, 2007.
  * Christakis and Fowler [2008] Nicholas A. Christakis and James H. Fowler.  The collective dynamics of smoking in a large social network.  _The New England Journal of Medicine_ , 358(21):2249–2258, 2008.  doi:10.1056/NEJMsa0706154.  URL https://doi.org/10.1056/NEJMsa0706154.
  * Supovitz et al. [2010] Jonathan Supovitz, Philip Sirinides, and Henry May.  How principals and peers influence teaching and learning.  _Educational administration quarterly_ , 46(1):31–56, 2010.
  * Rosenqvist [2018] Erik Rosenqvist.  Two functions of peer influence on upper-secondary education application behavior.  _Sociology of Education_ , 91(1):72–89, 2018\.
  * Sweet and Adhikari [2020] Tracy Sweet and Samrachana Adhikari.  A latent space network model for social influence.  _Psychometrika_ , 85(2):251–274, 2020.  doi:10.1007/s11336-020-09700-x.  URL https://doi.org/10.1007/s11336-020-09700-x.
  * Ogburn [2018] Elizabeth L Ogburn.  Challenges to estimating contagion effects from observational data.  _Complex spreading phenomena in social systems: influence and contagion in real-world social networks_ , pages 47–64, 2018.
  * McPherson et al. [2001] Miller McPherson, Lynn Smith-Lovin, and James M Cook.  Birds of a feather: Homophily in social networks.  _Annual review of sociology_ , 27(1):415–444, 2001.
  * Shalizi and Thomas [2011] Cosma Rohilla Shalizi and Andrew C. Thomas.  Homophily and contagion are generically confounded in observational social network studies.  _Sociological Methods & Research_, 40(2):211–239, 2011.  doi:10.1177/0049124111404820.  URL https://doi.org/10.1177/0049124111404820.  PMID: 22523436.
  * Spillane et al. [2018] James P. Spillane, Megan Hopkins, and Tracy M. Sweet.  School district educational infrastructure and change at scale: Teacher peer interactions and their beliefs about mathematics instruction.  _American Educational Research Journal_ , 55(3):532–571, 2018.  doi:10.3102/0002831217743928.  URL https://doi.org/10.3102/0002831217743928.
  * McFowland III and Shalizi [2021] Edward McFowland III and Cosma Rohilla Shalizi.  Estimating causal peer influence in homophilous social networks by inferring latent locations.  _Journal of the American Statistical Association_ , 0(0):1–12, 2021.  doi:10.1080/01621459.2021.1953506.  URL https://doi.org/10.1080/01621459.2021.1953506.
  * Cristali and Veitch [2022] Irina Cristali and Victor Veitch.  Using embeddings for causal estimation of peer influence in social networks.  In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, _Advances in Neural Information Processing Systems_ , volume 35, pages 15616–15628. Curran Associates, Inc., 2022.  URL https://proceedings.neurips.cc/paper_files/paper/2022/file/64587794695be22545d91c838243fcf8-Paper-Conference.pdf.
  * Liu et al. [2024] Jiewen Liu, Chan Park, Kendrick Li, and Eric J. Tchetgen Tchetgen.  Regression-based proximal causal inference, 2024.
  * Tchetgen [2014] Eric Tchetgen Tchetgen.  The control outcome calibration approach for causal inference with unobserved confounding.  _American Journal of Epidemiology_ , 179(5):633–640, 2014.  doi:10.1093/aje/kwt303.  URL https://doi.org/10.1093/aje/kwt303.
  * Miao et al. [2018] Wenbo Miao, Zhi Geng, and Eric Tchetgen Tchetgen.  Identifying causal effects with proxy variables of an unmeasured confounder.  _Biometrika_ , 105(4):987–993, 2018.  doi:10.1093/biomet/asy038.
  * Shi et al. [2020] Xi Shi, Wenyi Miao, and Eric Tchetgen Tchetgen.  A selective review of negative control methods in epidemiology.  _Current Epidemiology Reports_ , 7(4):190–202, 2020.  doi:10.1007/s40471-020-00243-4.  URL https://doi.org/10.1007/s40471-020-00243-4.
  * Wang Miao and Zhou [2023] Elizabeth L. Ogburn Wang Miao, Wenjie Hu and Xiao-Hua Zhou.  Identifying effects of multiple treatments in the presence of unmeasured confounding.  _Journal of the American Statistical Association_ , 118(543):1953–1967, 2023.  doi:10.1080/01621459.2021.2023551.  URL https://doi.org/10.1080/01621459.2021.2023551.
  * Louizos et al. [2017] Christos Louizos, Uri Shalit, Joris Mooij, David Sontag, Richard Zemel, and Max Welling.  Causal effect inference with deep latent-variable models.  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , NIPS’17, page 6449–6459, Red Hook, NY, USA, 2017. Curran Associates Inc.  ISBN 9781510860964.
  * Egami and Tchetgen Tchetgen [2023] Naoki Egami and Eric J Tchetgen Tchetgen.  Identification and estimation of causal peer effects using double negative controls for unmeasured network confounding.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 86(2):487–511, 12 2023.  ISSN 1369-7412.  doi:10.1093/jrsssb/qkad132.  URL https://doi.org/10.1093/jrsssb/qkad132.
  * Liu and Tchetgen Tchetgen [2022] Lan Liu and Eric Tchetgen Tchetgen.  Regression-based negative control of homophily in dyadic peer effect analysis.  _Biometrics_ , 78(2):668–678, 2022.  doi:https://doi.org/10.1111/biom.13483.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/biom.13483.
  * Dabbs et al. [2020] Beau Dabbs, Samrachana Adhikari, and Tracy Sweet.  Conditionally independent dyads (cid) network models: A latent variable approach to statistical social network analysis.  _Social Networks_ , 63:122–133, 2020.  ISSN 0378-8733.  doi:https://doi.org/10.1016/j.socnet.2020.06.004.  URL https://www.sciencedirect.com/science/article/pii/S0378873320300423.
  * Hoff et al. [2002] Peter D Hoff, Adrian E Raftery, and Mark S Handcock.  Latent space approaches to social network analysis.  _Journal of the American Statistical Association_ , 97(460):1090–1098, 2002.  doi:10.1198/016214502388618906.  URL https://doi.org/10.1198/016214502388618906.
  * Pitts and Spillane [2009] Virginia M. Pitts and James P. Spillane.  Using social network methods to study school leadership.  _International Journal of Research & Method in Education_, 32(2):185–207, 2009.  doi:10.1080/17437270902946660.  URL https://doi.org/10.1080/17437270902946660.
  * Forastiere et al. [2021] Laura Forastiere, Edoardo M. Airoldi, and Fabrizia Mealli.  Identification and estimation of treatment and interference effects in observational studies on networks.  _Journal of the American Statistical Association_ , 116(534):901–918, 2021.  doi:10.1080/01621459.2020.1768100.  URL https://doi.org/10.1080/01621459.2020.1768100.
  * Forastiere et al. [2022] Laura Forastiere, Fabrizia Mealli, Albert Wu, and Edoardo M. Airoldi.  Estimating causal effects under network interference with bayesian generalized propensity scores.  _Journal of Machine Learning Research_ , 23(289):1–61, 2022.  URL http://jmlr.org/papers/v23/18-711.html.
  * Aronow and Samii [2017] Peter M. Aronow and Cyrus Samii.  Estimating average causal effects under general interference, with application to a social network experiment.  _The Annals of Applied Statistics_ , 11(4):1912 – 1947, 2017.  doi:10.1214/16-AOAS1005.  URL https://doi.org/10.1214/16-AOAS1005.
  * Splawa-Neyman et al. [1990] Jerzy Splawa-Neyman, D. M. Dabrowska, and T. P. Speed.  On the Application of Probability Theory to Agricultural Experiments. Essay on Principles. Section 9.  _Statistical Science_ , 5(4):465 – 472, 1990\.  doi:10.1214/ss/1177012031.  URL https://doi.org/10.1214/ss/1177012031.
  * Rubin [1974] D.B. Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of Educational Psychology_ , 66(5):688–701, 1974.
  * Imbens and Rubin [2015] Guido W. Imbens and Donald B. Rubin.  _Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction_.  Cambridge University Press, USA, 2015.  ISBN 0521885884.
  * Fortunato [2010] Santo Fortunato.  Community detection in graphs.  _Physics Reports_ , 486(3):75–174, 2010.  ISSN 0370-1573.  doi:https://doi.org/10.1016/j.physrep.2009.11.002.  URL https://www.sciencedirect.com/science/article/pii/S0370157309002841.
  * Holland et al. [1983] Paul W. Holland, Kathryn Blackmond Laskey, and Samuel Leinhardt.  Stochastic blockmodels: First steps.  _Social Networks_ , 5(2):109–137, 1983.  ISSN 0378-8733.  doi:https://doi.org/10.1016/0378-8733(83)90021-7.  URL https://www.sciencedirect.com/science/article/pii/0378873383900217.
  * Karrer and Newman [2011] Brian Karrer and M. E. J. Newman.  Stochastic blockmodels and community structure in networks.  _Phys. Rev. E_ , 83:016107, Jan 2011.  doi:10.1103/PhysRevE.83.016107.  URL https://link.aps.org/doi/10.1103/PhysRevE.83.016107.
  * Hoff [2005] Peter D Hoff.  Bilinear mixed-effects models for dyadic data.  _Journal of the American Statistical Association_ , 100(469):286–295, 2005.  doi:10.1198/016214504000001015.  URL https://doi.org/10.1198/016214504000001015.
  * Chipman et al. [2010] Hugh A. Chipman, Edward I. George, and Robert E. McCulloch.  BART: Bayesian additive regression trees.  _The Annals of Applied Statistics_ , 4(1):266 – 298, 2010.
  * Lee and Song [2004] Sik-Yum Lee and Xin-Yuan Song.  Evaluation of the bayesian and maximum likelihood approaches in analyzing structural equation models with small sample sizes.  _Multivariate Behavioral Research_ , 39(4):653–686, 2004.  doi:10.1207/s15327906mbr3904_4.  URL https://doi.org/10.1207/s15327906mbr3904_4.  PMID: 26745462.
  * Sweet et al. [2013] Tracy M. Sweet, Andrew C. Thomas, and Brian W. Junker.  Hierarchical network models for education research: Hierarchical latent space models.  _Journal of Educational and Behavioral Statistics_ , 38(3):295–318, 2013.  doi:10.3102/1076998612458702.  URL https://doi.org/10.3102/1076998612458702.
  * Borg and Groenen [2005] Ingwer Borg and Patrick J. F. Groenen.  _Modern Multidimensional Scaling Theory and Applications_.  Springer, New York, 2005.  ISBN 038728981X.  doi:10.1007/0-387-28981-X.
  * Adhikari et al. [2021] Samrachana Adhikari, Tracy Sweet, and Brian Junker.  Analysis of Longitudinal Advice-Seeking Networks Following Implementation of High Stakes Testing.  _Journal of the Royal Statistical Society Series A: Statistics in Society_ , 184(4):1475–1500, 07 2021.  ISSN 0964-1998.  doi:10.1111/rssa.12708.  URL https://doi.org/10.1111/rssa.12708.
  * Spanbauer and Sparapani [2021] Charles Spanbauer and Rodney Sparapani.  Nonparametric machine learning for precision medicine with longitudinal clinical trials and bayesian additive regression trees with mixed models.  _Statistics in Medicine_ , 40(11):2665–2691, 2021\.  doi:https://doi-org.ezproxy.med.nyu.edu/10.1002/sim.8924.  URL https://onlinelibrary-wiley-com.ezproxy.med.nyu.edu/doi/abs/10.1002/sim.8924.
  * Hoshino [2008] Takahiro Hoshino.  A bayesian propensity score adjustment for latent variable modeling and mcmc algorithm.  _Computational Statistics & Data Analysis_, 52(3):1413–1429, 2008.  ISSN 0167-9473.  doi:https://doi.org/10.1016/j.csda.2007.03.024.  URL https://www.sciencedirect.com/science/article/pii/S0167947307001429.
  * Zigler et al. [2013] Corwin M. Zigler, Krista Watts, Robert W. Yeh, Yun Wang, Brent A. Coull, and Francesca Dominici.  Model feedback in bayesian propensity score estimation.  _Biometrics_ , 69(1):263–273, 2013.  doi:https://doi.org/10.1111/j.1541-0420.2012.01830.x.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1541-0420.2012.01830.x.
  * Kaplan and Chen [2012] David Kaplan and Jianshen Chen.  A two-step bayesian approach for propensity score analysis: Simulations and case study.  _Psychometrika_ , 77(3):581–609, 2012.  doi:10.1007/s11336-012-9262-8.
  * McCandless et al. [2010] Lawrence C McCandless, Ian J Douglas, Stephen J Evans, and Liam Smeeth.  Cutting feedback in bayesian regression adjustment for the propensity score.  _The International Journal of Biostatistics_ , 6(2), 2010\.  doi:10.2202/1557-4679.1205.
  * Bayarri et al. [2009] M. J. Bayarri, J. O. Berger, and F. Liu.  Modularization in Bayesian analysis, with emphasis on analysis of computer models.  _Bayesian Analysis_ , 4(1):119 – 150, 2009.  doi:10.1214/09-BA404.  URL https://doi.org/10.1214/09-BA404.
  * Warren et al. [2020] K. Warren, B. Campbell, and S. Cranmer.  Tightly bound: The relationship of network clustering coefficients and reincarceration at three therapeutic communities.  _Journal of Studies on Alcohol and Drugs_ , 81(5):673–680, September 2020.  doi:10.15288/jsad.2020.81.673.
  * Malik et al. [2016] N. Malik, F. Shi, H. W. Lee, and P. J. Mucha.  Transitivity reinforcement in the coevolving voter model.  _Chaos_ , 26(12):123112, December 2016.  doi:10.1063/1.4972116.
  * Hastie and Tibshirani [2000] Trevor Hastie and Robert Tibshirani.  Bayesian backfitting.  _Statistical Science_ , 15(3):196–213, 2000.  ISSN 08834237.  URL http://www.jstor.org/stable/2676659.

Figure 11: The teachers’ scores about student-centered beliefs for 9 schools
with more than 20 teachers in both 2012 and 2013. The belief scores are
represented visually by coloring each node based on high/low belief scores,
using a cut-off of $17$ for effective visualization.

##  Appendix A The detail of the real data

Item | Statment
---|---
1 | Teachers should encourage students to find their own solutions to math problems even if they are inefficient
2 | Teachers should allow students to figure out their own ways to solve simple word problems
3 | The goals of instruction in mathematics are best achieved when students find their own methods for solving problems
4 | Most students can figure out ways to solve many mathematics problems without any adult help
5 | Mathematics should be presented to children in such a way that they can discover relationships for themselves
Table 2: Teachers in Auburn Park were surveyed regarding their beliefs about
teaching mathematics. They were asked to indicate the extent to which they
agreed with the statements above.

##  Appendix B Proof for identification

Theorem 1 indicates that the average potential outcome $Y_{i2}(g)$ can be
obtained by an unbiased estimator of conditional outcome mean with proxy
variable $Z_{i}$. This identification enables us to estimate the causal peer
influence, which is determined by the function of potential outcome. Under
assumption 1-4,

$\displaystyle\operatorname{\mathbb{E}}\left[Y_{i2}(g)|\mathcal{A}_{n}\right]$
$\displaystyle=\operatorname{\mathbb{E}}\left[\operatorname{\mathbb{E}}\left[Y_{i2}(g)|Y_{i1}=y\operatorname{\bm{X}}_{i}=\operatorname{\bm{x}},U_{i}=u\right]|\mathcal{A}_{n}\right]$
$\displaystyle=\operatorname{\mathbb{E}}\left[\operatorname{\mathbb{E}}\left[Y_{i2}(g)|Y_{i1}=y,G_{i}=g,\operatorname{\bm{X}}_{i}=\operatorname{\bm{x}},U_{i}=u\right]|\mathcal{A}_{n}\right]\quad\text{by
Assumption \ref{assump:ignorability}}$
$\displaystyle=\operatorname{\mathbb{E}}\left[\operatorname{\mathbb{E}}\left[\operatorname{\mathbb{E}}\left[Y_{i2}(g)|Y_{i1}=y,G_{i}=g,\operatorname{\bm{X}}_{i}=\operatorname{\bm{x}},U_{i}=u\right]|Z_{i}=z\right]|\mathcal{A}_{n}\right]$
$\displaystyle=\operatorname{\mathbb{E}}\left[\operatorname{\mathbb{E}}\left[Y_{i2}(g)|Y_{i1}=y,G_{i}=g,\operatorname{\bm{X}}_{i}=\operatorname{\bm{x}},Z_{i}=z\right]|\mathcal{A}_{n}\right]\quad\text{by
tower property; $U_{i}$ is $Z_{i}$-measurable}$
$\displaystyle=\operatorname{\mathbb{E}}\left[\operatorname{\mathbb{E}}\left[Y_{i2}|Y_{i1}=y,G_{i}=g,\operatorname{\bm{X}}_{i}=\operatorname{\bm{x}},Z_{i}=z\right]|\mathcal{A}_{n}\right]\quad\text{by
Assumption \ref{assump:consistency} and Assumption
\ref{assump:Nei_interference}}$
$\displaystyle=\operatorname{\mathbb{E}}\left[Y_{i2}|Y_{i1}=y,G_{i}=g,\operatorname{\bm{X}}_{i}=\operatorname{\bm{x}},Z_{i}=z\right]\quad\text{by
Assumption \ref{assump:proxy}}$
$\displaystyle=\operatorname{\mathbb{E}}\left[Y_{i2}|Y_{i1}=y,G_{i}=g,\operatorname{\bm{X}}_{i}=\operatorname{\bm{x}},Z^{*}_{i}=z^{*}\right]\quad\text{by
$Z_{i}^{*}=WZ_{i}$ which is an invertible transformation}$

Given that $U$ is $Z$ measurable, tower property is
$E[X|Z]=E[E[X|Z]|U]=E[X|U]$. This proof adopts the proof presented by Cristali
and Veitch [2022], which also relies on the tower property.

##  Appendix C Prior specification and posterior computation in the semi-
parametric BART

For the tree structures, we use the default prior for the BART model. Our
Gibbs sampling algorithm uses a Metropolis-within-Gibbs strategy to
iteratively update ($\mathcal{T}_{i},\mathcal{M}_{t}$) for $t=1,\ldots,m$.
Sampling from ($\mathcal{T}_{i};\mathcal{M}_{t}$) is equivalent to sampling
from

$p\left(\mathcal{T}_{t},\mathcal{M}_{t}\mid\bm{R}_{t},\tau\right)$

where $\bm{R}_{t}=\left(R_{t1},\ldots,R_{tn}\right)^{\top}$ with
$R_{ti}=Y_{i}-\sum_{j\neq
t}g\left(D_{i}^{k};\mathcal{T}_{j},\mathcal{M}_{j}\right)-(\beta_{z}^{k})^{T}\bm{Z}_{i}^{k}$.
This allows us to use the existing Bayesian backfitting algorithm [Hastie and
Tibshirani, 2000, Chipman et al., 2010] to update the
$\left(\mathcal{T}_{t},\mathcal{M}_{t}\right)$.

We use multivariate Normal prior distribution $N_{d}(0,\sigma_{\beta}I)$ for
$\boldsymbol{\beta}^{k}_{z}$. The conditional posterior for the
$\boldsymbol{\beta}_{z}^{k}$ is

$\boldsymbol{\beta}_{z}^{k}\sim
N_{d}\left(\Psi^{k}\sum_{i=1}^{n_{k}}\bm{Z}_{i}^{k}\left(Y_{i2}^{k}-\mu_{i}^{k}\right)/\tau^{\prime
2},\Psi^{k}\right)$

where
$\Psi^{k}=\left(\sum_{i=1}^{n_{k}}\bm{Z}_{i}^{k}\bm{Z}_{i}^{kT}/\tau^{\prime
2}+(\sigma_{\beta}I)^{-1}\right)^{-1}$ and $\mu_{i}^{k}=Y_{i}-\sum_{j\neq
t}g\left(D_{i}^{k};\mathcal{T}_{j},\mathcal{M}_{j}\right)$.

Algorithm 2 A single iteration for HLSM-BART

1: (Latent location stage)

2:for $k=1$ to $K$ do

3: Draw a sample of HLSM parameters $\theta_{z}^{k}\sim
p(\cdot|\mathbf{Z}^{k},A^{k})$ where
$\mathbf{Z}^{k}=(Z_{1}^{k},\cdots,Z_{n}^{k})$

4: for $i=1$ to $n_{k}$ do

5: Update $Z_{i}^{k}$ given $Z_{i}^{k},\theta_{z}^{k}$ and $A$

6: end for

7:end for

8: (Outcome stage)

9:for $m=1$ to $M$ do

10: Set $R^{k}_{mi}=y^{k}_{i}-\sum_{q\neq
m}g\left(\bm{x}^{k}_{i};\mathcal{T}_{q},\mathcal{M}_{q}\right)-(\beta_{z}^{k^{\prime}})^{T}\bm{Z}^{k}_{i}$
for $i=1,\ldots,n$ and $k=1,\ldots,K$.

11: Propose a new tree $T_{m}^{*}\sim
Q\left(\mathcal{T}_{m}\rightarrow\mathcal{T}_{m}^{*}\right)$.

12: Compute the acceptance ratio  $\displaystyle
r=\frac{\mathscr{L}\left(\mathcal{T}_{m}^{*}\right)Q\left(\mathcal{T}_{m}^{*}\rightarrow\mathcal{T}_{m}\right)}{\mathscr{L}\left(\mathcal{T}_{m}\right)Q\left(\mathcal{T}\rightarrow\mathcal{T}_{m}^{*}\right)},$
where
$\mathscr{L}\left(\mathcal{T}_{m}\right)=p\left(\mathcal{T}_{m}\mid\bm{R}_{m},\sigma\right)\times
p\left(\mathcal{T}_{m}\right)$ and
$\bm{R}_{m}=(R^{1}_{m1},\ldots,R^{1}_{mn_{1}},\ldots,R^{k}_{m1},\ldots,R^{k}_{mn_{k}})$.

13: Set $\mathcal{T}_{m}=\mathcal{T}_{m}^{*}$ if $U\leq\min(1,r)$ where
$U\sim\operatorname{Uniform}(0,1)$, otherwise retain $\mathcal{T}_{m}$.

14: Sample the $\mu_{t\ell}$’s from their full conditional distributions.

15:end for

16:Sample $\tau^{\prime}$ from its full conditional distribution.

17:Sample $\beta_{z}^{k^{\prime}}$ from its full conditional distribution.

18: (Monte Carlo estimation stage)

19:for $g\in\mathcal{G}$ do

20: Predict $\hat{Y}^{k}_{i2}(g)$ given $Z^{k}_{i}$,
$\bm{\mathcal{T}}=(\mathcal{T}_{1},\ldots,\mathcal{T}_{m}),$ and
$\bm{\mathcal{M}}=(\mathcal{M}_{1},\ldots,\mathcal{M}_{m})$.

21:end for

22:Average the potential outcomes over all units
$\hat{\mu}(y,g)=\frac{1}{NK}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}\hat{Y}^{k}_{i2}(g)$

##  Appendix D Details of data generating process and estimators for the
multiple network setting

Given the generated multiple networks following HLSM in (10), the individual
behaviors at baseline are assigned as following

$\displaystyle
Y_{i1}^{k}=5+\beta_{u_{1}}U_{i1}^{k}+\beta_{u_{2}}U_{i2}^{k}+\epsilon_{i1}^{k}$
(15)

where $\epsilon_{i1}^{k}\overset{iid}{\sim}N(0,1)$. We further specify a
nonlinear outcome model for individual behaviors at follow-up as

$\displaystyle
Y_{i2}^{k}=10+3Y_{i1}^{k}+10\text{sigmoid}(10(G_{i}^{k}-5))+\beta^{\prime}_{u_{1}}U_{i1}^{k}+\beta^{\prime}_{u_{2}}U_{i2}^{k}+\epsilon_{i2}^{k}$
(16)

where $\text{sigmoid}(x)=1/(1+e^{-x})$ and
$\epsilon_{i2}^{k}\overset{iid}{\sim}N(0,1)$. Note that
$G_{i}^{k}=\sum_{j\in\mathcal{N}^{k}_{i}}Y_{j1}^{k}/\sum_{j\in\mathcal{N}^{k}_{i}}A_{ij}$
where $\mathcal{N}^{k}_{i}$ is a set of neighbors that share a tie with node
$i$ in the $k$th network. The parameter values are fixed at
$\beta_{\mathbf{u}}=(\beta_{u_{1}},\beta_{u_{2}})=(0.5,0.5)$ and
$\beta^{\prime}_{\mathbf{u}}=(\beta^{\prime}_{u_{1}},\beta^{\prime}_{u_{2}})=(10,10)$.
With these true underlying values, the population average potential outcome is
$\operatorname{\mathbb{E}}(Y_{i2}(g))=25+10\text{sigmoid}(10(g-5))$.

In the estimation step to estimate latent locations, we set hyperparameters as
$\sigma^{\prime}_{z}=5,\ \mu_{\alpha^{\prime}}=0,\ \mu_{\beta^{\prime}}=0,\
\sigma_{\alpha^{\prime}}^{2}=100,\ \sigma_{\beta^{\prime}}^{2}=100$ in (10).
Similar to the single network setting, three different estimators are
compared, each using a distinct set of predictors in the BART model:

  * •

Oracle: $Y^{k}_{i2}=f(D_{i}^{k})+\epsilon_{i}^{k}$ where
$\mathcal{D}^{k}_{i}=(Y^{k}_{i1},G^{k}_{i},U^{k}_{i1},U^{k}_{i2})$.

  * •

LSM-semiBART:
$Y^{k}_{i2}=f(D_{i}^{k})+\beta_{z_{1}}^{k}Z_{i1}^{k}+\beta_{z_{2}}^{k}Z_{i2}^{k}+\epsilon_{i}^{k}$
where $\mathcal{D}^{k}_{i}=(Y^{k}_{i1},G^{k}_{i})$.

  * •

Naive: $Y^{k}_{i2}=f(D_{i}^{k})+\epsilon_{i}^{k}$ where
$\mathcal{D}^{k}_{i}=(Y^{k}_{i1},G^{k}_{i})$.

where $Z_{i1}^{k}$ and $Z_{i2}^{k}$ are estimated latent locations from HLSM
in (10) and $f(\cdot)$ is modeled as the BART model.

##  Appendix E Supplementary plots for simulation and real data analysis

Figure 12: Simulation 1: The generated networks from the data generating
process with $\sigma_{j}\in(0.1,\ldots,1)$ (top panels) and corresponding
posterior mean of latent locations (bottom panels) Figure 13: Simulation 3:
The latent positions, $U_{i}$’s, from data generating process (left) and
posterior mean of latent positions $Z_{i}$’s (right) corresponding to $1^{st}$
and $5^{th}$ network. Figure 14: The estimated average belief scores in 2013
across all school from the Naive (left) and LSM-semiBART (right) estimator is
represented for each school.

##  Appendix F goodness of fit of HLSM

In this section, we evaluate the goodness of fit of HLSM using advice-seeking
network data. We perform posterior predictive checks by generating 11 networks
from the parameters obtained at each MCMC iteration. Specifically, we utilize
the latent space positions and the random intercepts $\alpha_{0}^{k}$ to
generate these networks. We assess the model fit based on network density and
transitivity, which can partially explain homophily. We use 200,000 MCMC
iteration with 5000 burn-in thinning by every 200th step.

First, we visualize an example of a replicated network alongside the observed
network in Figure 15. The ties of both the observed and replicated networks
are depicted on the estimated latent positions from a randomly selected MCMC
iteration. Generally, this visualization demonstrates that the network ties
from the replicated networks closely resemble those of the observed network.

Figure 16 illustrates the distributions of network (tie) density and
transitivity across the generated networks. The observed network statistic for
each network are shown with a vertical line. In general, density and
transitivity are adequately recovered indicating that the HLSM was able to
capture these characteristics of the observed network. In School 2, however,
both density and transitivity are overestimated potentially due to a high
average degree (3.6), which results in closely positioned latent locations,
thereby overestimating density and transitivity. To address this limitation
and more effectively capture complex network characteristics, considering
different dimensionalities of $Z$ for each school could be explored for future
work.

Figure 15: The illustration of observed and replicated network ties on
estimated latent positions from a randomly selected MCMC sample. Figure 16:
Posterior predictive checks of the HLSM: tie density (upper panel) and
transitivity (lower panel). Figure 17: Density plot of the likelihood values
computed at each MCMC iteration (left) and the posterior distributions of
$\alpha_{0}^{1}$ (right) derived from advice-seeking networks for each
thinning (1,20,100,200). Each plot includes the posterior mean and the 95%
credible intervals, depicted on the bottom line and annotated within the plot.