Joint mixed-effects models for causal inference in clustered network-based observational studies

1. Introduction

Much of the theoretical development in causal inference has been within the counterfactual framework, which typically employs a ‘no interference’ assumption, in that the potential outcome of an individual is not affected by the treatments of another. In many sociological settings, the ‘no interference’ assumption is implausible, and traditional methods will fail whenever social interaction influences the mechanisms relating a treatment to an outcome. In such cases, the treatment or exposure assigned to an individual will also affect their close neighbours’ outcomes, yielding spillover or indirect effects. A highly relevant example is in studies of adolescents, where the influence of peers takes on a much greater role than in early childhood or later life. For instance, access to drugs in the home (exposure) may affect that teen’s HIV risk and that of their friends, as adolescents are particularly susceptible to peer influence. A tobacco cessation intervention in a high school setting might not only affect the uptake of its direct recipients but also that of their contacts through word of mouth and social influence. Similarly, it is entirely plausible that an adolescent’s academic performance is influenced both by their home environment and that of their peers.

Peer influence has long been of interest to social scientists, since it can potentially inform the optimal organisation of schools, neighbourhoods, workplaces, and other structures whereby individuals interact. In educational research, there has been mounting evidence suggesting that classroom composition, notably in terms of ability or achievement distribution, can have an influence on learning outcomes.^1,2 For instance, using administrative public school data gathered by the Texas Education Agency, Hoxby ³ found that a change in one point in peers’ reading scores is associated with an increase in a student’s own score between 0.15 and 0.40. Leveraging data from the National Longitudinal Study of Adolescent Health (Add Health), Bifulco et al. ⁴ found that having a higher proportion of classmates with college-educated mothers was associated with a lower risk of dropping out and a higher likelihood of attending college. Bifulco et al. ⁴ also hypothesised about a spillover effect of maternal educational attainment on students’ post-secondary health behaviors, such as binge drinking, cigarette smoking, and marijuana use, but the data did not provide evidence of such an effect. Fletcher et al. ⁵ made use of the friendship nominations data in Add Health and found that a greater proportion of friends with college-educated mothers is associated with a higher grade point average (GPA), with evidence of a differential effect by sex. This body of evidence supports the idea of spillover effects of home environment on academic performance, and possibly health behaviour outcomes.

Social network data are often collected when such interference is at play, and many of the settings in which spillover effects are being estimated are multilevel (e.g. social networks of students nested within schools). The Add Health study is one example of such a multilevel structure. When faced with multilevel network data, such as one might find in schools or geographical communities, for example, there may be an unmeasured, latent factor influencing both the exposure and the outcome at the cluster level. This tendency for peers, such as schoolmates or neighbourhood residents, to share some unobserved contextual factors is termed contextual confounding. ⁶ In the context of adolescent school performance, there could be a latent trait influencing both the average maternal education as well as the adolescents’ overall performance at the school level. If it is unmeasured, such a latent trait could be geographical location, for instance, since it may encompass area-level violence, pollution, or access to quiet learning spaces. Another potential source of network confounding is homophily, which is defined as the increased tendency of units with similar characteristics to form ties. As opposed to contextual confounding, homophily confounding arises at the individual rather than the cluster level and pertains to the network structure within clusters. In this paper, the main focus is on unmeasured contextual confounding, and latent homophily will either be assumed to be not present or removed by conditioning on all the variables driving tie formation in simulations.

Causal inference is inherently a missing data problem, and the problem of unmeasured contextual confounding can be conceived as a non-ignorable missingness mechanism. ⁷ Non-ignorable missingness or treatment assignment mechanisms are typically addressed by imposing additional parametric assumptions. In the frequentist paradigm, joint mixed-effects models (JMMs) were originally proposed to handle non-ignorable missing data when estimating the parameters of an outcome data-generating model.^8,9 Shardell and Ferrucci ¹⁰ combined parametric JMMs for a study outcome and an exposure with g-computation to identify and estimate causal effects in longitudinal settings under unmeasured confounding. In the Bayesian paradigm, Papadogeorgou and Samanta ¹¹ proposed an approach based on simultaneous modelling of the exposure and outcome to account for the presence of spatially structured unmeasured confounders. Nobre et al. ¹² used a Bayesian multilevel approach to study the effect of directly observed therapy on the rate of cure of tuberculosis while mitigating the impact of unobserved contextual factors shared by individuals within different municipalities across the state of São Paulo, Brazil. Xu et al. ¹³ proposed Bayesian multivariate generalised linear mixed-effects models (LMMs) to assess dynamic treatment regimes while accounting for unmeasured patient heterogeneity. The aforementioned methods employ an outcome modelling strategy while leveraging the treatment assignment to learn about the cluster-level confounders. Other methods based on propensity score modelling have been proposed to reduce bias from unmeasured cluster-level confounding,^14–17 but these approaches are suitable for multilevel designs where there are no spillover effects through interference.

Existing methods that accommodate multilevel networks, be it inverse probability-of-treatment weighting, outcome regression, or doubly robust estimation, predominantly rely upon the assumption of no unmeasured confounding.^18–20 In this work, we address an important gap by relaxing the unconfoundedness assumption and providing a Bayesian framework for causal inference in the multilevel network setting, in which heterogeneity is due to clustering of individuals in different subgraphs. We propose an approach based on a JMM combined with Bayesian standardisation to estimate direct and indirect effects of a binary treatment on a continuous outcome in the presence of unmeasured contextual confounding. We study frequentist properties of the JMM approach in simulation studies, which show that valid inference can be obtained in comparison with more naive approaches that either assume ignorability of the treatment assignment (LMM) or ignore the clustered nature of the data altogether (fixed effects model (FEM)).

The paper is structured as follows. In Section 2, we introduce the Add Health study and our motivating question examining the direct and indirect effect of maternal education on school performance. In Sections 3 and 4, we introduce the counterfactual framework and causal assumptions under the network interference setting and the relevant causal estimands for this setting. In Section 5, we describe the proposed JMM approach along with the Bayesian standardisation procedure used to perform posterior inference for causal effects. In Section 6, we assess frequentist properties of the proposed estimator in simulation studies. The method is then utilised to assess the direct and spillover effects of maternal college education on adolescent school performance using the Add Health network data in Section 7. We close with a general discussion in Section 8.

2. The Add Health study

Add Health is a longitudinal study of a nationally representative sample of 90,118 adolescent students in grades 7 to 12 in the United States in 1994–1995 who were followed throughout adolescence. ²¹ The original purpose of the study was to understand the determinants of adolescent health and health behaviour. Add Health used a school-based design that selected 80 high schools and associated feeder schools (typically a middle school) using a sampling frame derived from the Quality Education Database. The final sample comprised 132 schools with communities located in urban, suburban, and rural areas of the country.

In addition to socio-demographic and academic achievement characteristics, peer network data were obtained in the in-school questionnaire, in which adolescents nominated up to five best male and five best female friends from the school roster. We mapped an undirected network to friendship nomination data from the first wave of data collection in 1994–1995 of the Add Health. Following a reciprocal definition of friendship, we drew an edge between students i and j if student i listed j as a friend in the in-school survey, or j listed i as a friend, or both. In other words, friendships that were not mutual were retained in the graph. For the purpose of this analysis, we selected a random sample of 16 schools, which allowed us to reduce the computational burden. Figure 1 represents the three smallest schools in the subsample, with students with a college-educated mother represented by black nodes. The network characteristics and distribution of attributes are summarised in Tables S4 and S5, respectively, which can be found in Supplemental Appendix D. The Add Health subnetwork has an average degree of 6.65 (standard deviation = 3.88) and a density of 0.0008. The transitivity, which measures the extent to which nodes tend to cluster together, is 0.171. Perfect transitivity implies that if two nodes x and z are connected to another node y , then systematically x and z are also connected. The assortativity by maternal education, which is 0.123 globally, is the extent to which nodes with matching values of maternal education tend to form ties in the network.

OPEN IN VIEWER

Figure 1.

Three schools in the Add Health network of 71, 159, and 160 students, respectively, where an edge between two students indicates that either one nominated the other as a friend, and black nodes represent students with college-educated mothers. Note that ties that exist across different schools are not shown in this visualisation.

Previous work has examined the causal effect of peers on academic success under a counterfactual framework, focusing on the indirect (spillover) effects of having a college-educated mother on GPA. ²⁰ As individual-level data are available for each school in Add Health, it is expected that students attend the same school and thus reside in the same area share latent traits affecting both their propensity to have a college-educated mother and their overall grades. Thus, we suspect unmeasured subgraph-level confounders influencing both the overall socio-economic status and the overall academic performance at the school level. This raises the question of whether subgraph-level random effects should be included when modelling the quantities of interest. Therefore, we follow the Bayesian paradigm and discuss a procedure that jointly models propensity score and the outcome to estimate direct and indirect effects of maternal education on academic performance. Throughout the text, we will refer to maternal education as a (self-selected) treatment.

Consider an undirected network G = ( N , E ) , where N = { 1 , 2 , … , N } is a set of N nodes and E is a set of edges E with generic element { i , j } = { j , i } denoting the presence of an edge between nodes i and j . The sociomatrix A associated with G has an entry A i j = 1 if { i , j } is in E , and 0 otherwise. We assume that the set of nodes in this network is clustered in that it can be written as a partition of m sets, such that N = ∪ ν = 1 m N ν , where N ν ∩ N υ = ∅ for ν ≠ υ and each set is of order | N ν | = N ν . The m disjoint node sets could, for instance, correspond to distinct schools or classrooms. We define G ν = ( N ν , E ν ) to be the induced subgraph associated with the set of nodes N ν . Note that this definition does not preclude the existence of ties across the different disjoint node sets N ν , which implies that the induced subgraphs { G ν , ν = 1 , … , m } are not necessarily disjoint. We define the neighbourhood N ν i of node i in cluster ν , i = 1 , 2 , … , N ν , as the set of nodes j for which { i , j } ∈ E ν and the degree of node i is denoted d ν i = | N ν i | . Since the induced subgraphs { G ν , ν = 1 , … , m } are not necessarily disjoint, nodes in the immediate neighbourhood N ν i may belong to a different subset than N ν , say N υ , υ ≠ ν . Additionally, we denote N ν i * = N ν i ∪ i as the union of the node i and the nodes in N ν i . It follows that each node i is associated with a partition of N , that is, { i , N ν i , N ν − i } , where N ν − i = N ∖ N ν i * .

A node i in cluster or subgraph ν has an observed outcome Y ν i , an assigned binary treatment Z ν i , and a vector of covariates X ν i . Covariates may include an intercept, unit i ’s individual covariates, or covariates corresponding to unit i ’s neighbourhood or subgraph. Let Z N ν i = ( Z j 1 , … , Z j d ν i ) ′ denote the observed neighbourhood treatment vector for node i in subgraph ν , where j k ∈ N ν i , k = 1 , … , d ν i . We extend typical causal assumptions to include the presence of network interference. Let y ν i ( z ν i , z N ν i , z N ν − i ) denote the potential outcome for individual i in subgraph ν if they (and their neighbours and non-neighbours) receive the treatment assignment ( z ν i , z N ν i , z N ν − i ) . While in principle this allows the potential outcome to depend on the entire treatment assignment in the network, we will soon restrict the number of potential outcomes using the known individual connections and assuming an exposure mapping; see Assumption 2 below.^19,22 We first make the consistency assumption, which states that there cannot be multiple versions of a treatment.

Assumption 2 (Neighbourhood interference)

For ν ∈ { 1 , … , m } , i ∈ N ν , ∀ z N ν − i , z N ν − i ′ , and ∀ z N ν i , z N ν i ′ such that ϕ ( z N ν i ) = ϕ ( z N ν i ′ ) given a function ϕ : { 0 , 1 } d ν i → Φ ν i , then,

y ν i ( z ν i , z N ν i , z N ν − i ) = y ν i ( z ν i , z N ν i , z N ν − i ′ )

This assumption is akin to Assumption 2 from Forastiere et al., ¹⁹ which is that the outcome for an individual depends on the exposures of other individuals only through a function of the neighbourhood treatment vector. For simplicity, in what follows we assume that the function ϕ corresponds to the total exposure value among the neighbours, that is, ϕ : { 0 , 1 } d ν i → Φ ν i defined by ϕ ( z N ν i ) = ∑ j ∈ N ν i z ν j , where Φ ν i ∈ { 0 , 1 , … , d ν i } . Because of Assumption 2, we can drop the dependence of y ν i ( ⋅ ) on the nodes in N ν − i and write more simply y ν i ( z ν i , z N ν i ) .

In the presence of unmeasured contextual confounding, the covariates X ν i are no longer a sufficient conditioning set for unconfoundedness of the treatment assignment, that is, y ν i ( z ν i , z N ν i ) Z ν i , Z N ν i | X ν i . We assume, however, that exchangeability holds conditional on measured covariates and the subgraph-level value of the unmeasured covariate. For a given subgraph ν , we represent these subgraph-level unmeasured covariates as b ν = ( b ν y , b ν z ) , where b ν y is an unmeasured covariate pertaining to the outcome data generating mechanism, while b ν z informs the treatment assignment mechanism. The following assumption states that, for node i in subgraph ν , potential outcomes and treatment assignments are independent conditional on the value of the measured covariate vector, X ν i , and the unmeasured confounder b ν y .

Assumption 3 (Latent ignorability)

There exists unmeasured subgraph-level covariates b y = ( b 1 y , b 2 y , … , b m y ) , each of dimension q 1 , such that for all ν ∈ { 1 , … , m } , i ∈ N ν , z ν i ∈ { 0 , 1 } , z N ν i ∈ { 0 , 1 } d ν i ,

y i ( z ν i , z N ν i ) ⊥ ⊥ Z ν i , Z N ν i | X ν i , b ν y

Although inference in this work relies primarily upon the outcome model, we must also ensure appropriate covariate overlap, which tends to reduce the sensitivity of causal estimates to correct model specification ²³ and further ensures that the estimated treatment effect is not extrapolated to subregions of the covariate space for which there is little or no information. As the propensity score is central to the overlap or positivity assumption, we define the joint propensity score as e ( z ν i , z N ν i | X ν i , X N ν i , b ν z ) = P ( Z ν i = z ν i , Z N ν i = z N ν i | X ν i , X N ν i , b ν z ) , which corresponds to the probability of observing the joint treatment ( z ν i , z N ν i ) conditional on measured individual and neighbourhood covariates and the unmeasured subgraph-level covariate b ν z . Thus, we further require:

Assumption 4 (Positivity)

For ν ∈ { 1 , … , m } , i ∈ N and for all z ν i , X ν i , X N ν i , b ν z ,

e ( z ν i , z N ν i | X ν i , X N ν i , b ν z ) > 0

4. Estimands

To define the causal estimands of interest, we adopt the Bernoulli exposure allocation standardisation.^18,24 Let π ( Σ z N ν i ; α ) = ( d ν i Σ z N ν i ) α Σ z N ν i ( 1 − α ) d ν i − Σ z N ν i denote the probability of node i in subgraph ν receiving neighbourhood treatment Σ z N i under Bernoulli allocation strategy or treatment coverage α , which corresponds to the counterfactual probability of directly receiving treatment. In this counterfactual scenario, potential outcomes are standardised with respect to a population in which all individuals in a neighbourhood are exposed to treatment independently with the same probability α , which corresponds to a type B parameterisation. ²⁵ While independent exposure selection among neighbours seems plausible in the context of maternal education in Add Health, note that generalisations of the type B estimands have been proposed to allow for within-cluster dependence in treatment uptake in contexts where it would be warranted, such as those of vaccination and air pollution.^26,27 Define the i -th individual’s average potential outcome given individual exposure z and allocation strategy α conditional on covariate value x as

μ z α ν i ( x ) = E [ y ¯ ν i ( z ; α ) | X ν i = x ] = E [ ∑ Σ z N ν i = 0 d ν i y ν i ( z , Σ z N ν i ) π ( Σ z N ν i ; α ) | X ν i = x ]

We consider various contrasts of the average potential outcomes to define the causal estimands.^18,20,24,28 The direct exposure effect is defined as the difference between the average potential outcomes of untreated and treated individuals for a fixed allocation strategy α , that is, DE ( α ) = μ 1 α − μ 0 α . The spillover effect, also referred to as the indirect effect, is defined as the difference between the average potential outcomes under counterfactual scenarios α and α ′ among the untreated, that is, IE ( α , α ′ ) = μ 0 α − μ 0 α ′ . Note that it would be possible to define a spillover effect for the treated individuals as well. For allocation strategies α and α ′ , we also define the total effect TE ( α , α ′ ) = μ 1 α − μ 0 α ′ as well as the overall effect OE ( α , α ′ ) = μ α − μ α ′ . As the name suggests, TE ( α , α ′ ) can be interpreted as the total effect of treatment, since TE ( α , α ′ ) = DE ( α ) + IE ( α , α ′ ) , whereas the overall effect corresponds to the difference in average potential outcomes under one allocation strategy compared to another allocation strategy.

A summary of the different causal estimands presented along with the distributions being averaged over can be found in Table S1 in Supplemental Appendix A.

5. Proposed JMM

The proposed joint model allows for the possibility of unmeasured contextual confounding by assuming that the latent factors influencing the study outcome and the exposure selection in the social network are correlated, as illustrated in the directed acyclic graph (DAG) in Figure 2.

OPEN IN VIEWER

Figure 2.

Directed acyclic graph (DAG) for the joint linear mixed-effects model displaying interference and unmeasured heterogeneity in both treatment assignment and outcome mechanisms within a subgraph. Baseline characteristics X ν i and X ν j are excluded from the figure for simplicity.

5.1. Outcome model

We assume that the potential outcomes arise according to

y ν i ( z ν i , Σ z N ν i ) = β Z z ν i + h 1 ( Σ z N ν i ; β Z N ) + h 2 ( z ν i Σ z N ν i ; β Z × Z N ) + h 3 ( X ν i y ; β X ) + W ν i ⊤ b ν y + δ ν i

(1)

where h 1 , h 2 are known functions, X ν i y is a vector of nodal covariates that are contained in X ν i , W ν i is a vector of random effect covariates, and δ ν i is the random error of node i in subgraph ν . We further assume that δ = ( δ 1 , δ 2 , … , δ m ) ∼ N ( 0 , Φ ) , where δ ν = ( δ ν 1 , … , δ ν N ν ) , ν = 1 , 2 , … , m . In related work, the function h 1 is often assumed to be linear,^18,19 though non-linear functions such as spline representations could easily be accommodated. Similarly, h 2 and h 3 could either be linear or non-linear. According to this model, β Z and (potentially vector-valued) β Z N are tied to the direct and spillover effects of the exposure, respectively. The parameter β Z × Z N captures the interaction between z and Σ z N ν i and modifies the direct and interference effects with the level at which we fix the individual or neighbourhood exposure. In the absence of interactions between the exposure variables and the covariates/random effects, the causal effects do not vary as a function of characteristics of the node or the subgraph.

We assume a simultaneous autoregressive (SAR) model for the vector of random errors δ , such that if U = { I N − τ A } − 1 exists, the covariance matrix of δ is given by Cressie ²⁹

Φ = σ ε 2 U U ⊤

(2)

Due to its dependence on the adjacency matrix A , the SAR covariance matrix is such that the potential outcomes of neighbouring individuals are more correlated than those of non-neighbours. ³⁰ According to this model, the parameter τ controls the strength of the interaction between neighbours. ³¹ A simplifying assumption we make in simulations is that the subgraphs G ν are disjoint, such that the sociomatrix A associated with G is a block-diagonal matrix with m blocks denoted A ν , ν = 1 , … m , each corresponding to a subgraph G ν . This implies that Φ is also block-diagonal with blocks Φ ν = σ ε 2 U ν U ν ⊤ , ν = 1 , … , m , assuming that U ν = { I N ν − τ A ν } − 1 exists for all ν . If this assumption is not met, the computation might be infeasible in the presence of a large graph order N due to difficulties with large-scale matrix inversion.

5.2. Exposure model

We posit a Bernoulli distribution for the individual assigned treatment, such that Z ν i | p ν i ∼ Bernoulli ( p ν i ) with

logit ( p ν i ) = f ( X ν i z ; γ ) + U ν i ⊤ b ν z

(3)

where p ν i = P ( Z ν i = 1 | X ν i z , b ν z ) , f is a known function, and U ν i is a vector of random effect covariates. Note that the vector of covariates X ν i z in the exposure selection model need not be identical to X ν i y . In a standard causal setting, it is typically preferable to avoid including instruments in the treatment model, notably when implementing inverse weighting approaches, since it can lead to extreme weights and decreased precision in the estimation of causal effects. ³² However, since we do not assume conditional exchangeability, a correctly specified model that follows the data-generating procedure (accounting for more than simply the confounders) for the treatment assignment mechanism is required for unbiased inference. Thus, X ν i z is not necessarily a subset of X ν i y in the proposed joint modelling methodology.

In what follows, similar to Lee et al. ¹⁸ and McNealis et al., ²⁰ we make the simplifying assumption that conditional on the cluster-level random effect b ν z , for ν = 1 , … , m , i = 1 , … , N ν , the joint propensity score can be expressed as

e ( z ν i , z N ν i | X ν i , X N ν i , b ν z ) = P ( Z ν i = z ν i , Z N ν i = z N ν i | X ν i , X N ν i , b ν z ) = ∏ j ∈ { i , N i } P ( Z ν j = z ν j | X ν j , b ν z )

When describing the full data likelihood in Section 5.4, this assumption allows us to express the likelihood for the exposure model parameters as the product of individual propensity scores.

5.3. Random effect distribution

The subgraph-level random effects of the above models, b ν = ( b ν y , b ν z ) , are a q 1 + q 2 dimensional vector that captures the unobserved dependence among observations within each subgraph. We denote b = ( b 1 , … , b m ) the collection of subgraph-level random effects across the whole graph G . We assume that the random effects follow a multivariate normal distribution b ν ∼ N q 1 + q 2 ( 0 , Σ ) , where

Σ = ( Σ 11 Σ 12 Σ 21 Σ 22 )

(4)

The correlations of the random effects in the two models represent unobserved dependence among the outcome and the exposure assignment mechanisms. Groups defined by different treatment assignments are not, in the presence of such latent effects, conditionally exchangeable given observed covariates. This means that merely incorporating random effects in the outcome regression and ignoring the treatment assignment will not lead to valid inference of causal effects. Although based on Figure 2 and Assumption 3, conditioning on b y would appear to block confounding, these are unobserved covariates and their posterior distribution may only be accurately estimated by accounting for the dependence between the outcomes and treatment assignments, the latter being informative for the imputation of missing potential outcomes. ⁸ In other words, latent ignorability conditional on b y can only be achieved if b y is modelled in conjunction with b z in a joint inference procedure which accounts for the dependence in the outcome and treatment-specific random effects. A model that estimates the random effects in the outcome without accounting for their correlation with the random effects of the treatment model will be poorly estimated, and this measurement error in the latent variables may lead to bias in the treatment effect estimators.

The use of a latent Gaussian structure to relate the outcome and exposure assignment mechanisms is similar to the joint modelling strategies used by Tseng et al. ³³ in the context of non-ignorable missing data and Xu et al. ¹³ in the context of non-ignorable exposure assignment in longitudinal data. A different strategy, which is arguably more restrictive and more challenging to fit in comparison, is the use of a shared parameter model as the working model, where the outcome and the model for the exposure share common random effects. Examples of their use include Follmann and Wu ⁸ in the context of non-ignorable missingness and Shardell and Ferrucci ¹⁰ in the context of non-ignorable treatment assignment.

5.4. Bayesian inference

With interference, only one of multiple possible potential outcomes is observed for each node, with the number of potential outcomes for one unit depending on its degree by virtue of Assumption 2. Under the Bayesian paradigm, unobserved potential outcomes are viewed as missing data, and inference about causal estimands is obtained from their posterior distributions.^7,11,23,34 Let y ( ⋅ ) = { y i ( ⋅ ) for all i ∈ N } denote the collection of potential outcomes for all nodes in G and possible individual and neighbourhood exposures. We can then define the partition y ( ⋅ ) = ( Y , y ~ ) , where Y = ( Y 1 , Y 2 , … , Y m ) is the collection of observed outcomes across all the subgraphs and y ~ is the collection of missing potential outcomes for all graph units. We denote the collection of covariate vectors for all the nodes in G as X . As we do not assign a probability distribution for the multi-dimensional X , we consider a conditional likelihood in what follows. Following the specifications discussed in the previous section, the joint distribution of ( y ( ⋅ ) , Z ) conditional on X , denoted p ( y ( ⋅ ) , Z | X , θ ) , is governed by the parameter vector θ = ( β , σ ε 2 , τ , γ , Σ ) . Further, we can distinguish the parameters of the outcome model, θ Y = ( β , σ ε 2 , τ , Σ ) , and those of the exposure model, θ Z = ( γ , Σ ) , where Σ is shared between the two models. Bayesian causal inference proceeds by specifying a prior p ( θ ) and imputing missing potential outcomes from the posterior predictive distribution

p ( y ~ | O ) = ∫ θ p ( y ~ | θ , O ) p ( θ | O ) d θ ∝ ∫ θ [ ∫ b ν p ( y ~ | b , θ , O ) p ( b ν | θ , O ) d b ] L ( θ | O ) p ( θ ) d θ

where p ( y ~ | θ , O ) is the sampling distribution of the missing potential outcomes given the parameters and the observed data, p ( θ | O ) is the posterior distribution of θ , L ( θ | O ) = p ( O | θ ) is the observed data likelihood, and O = ( Y , Z , X ) is the observed data. Since our interest lies in marginal causal estimands, we perform inference from the marginal posterior predictive distribution of the potential outcome with respect to the random effects, which is highlighted by the proportionality statement.

Typically, the central component in Bayesian inference for causal effects is the specification of the outcome model. ²³ In instances where the propensity score model is ignorable, unbiased inference can be obtained about causal inference without having to specify a model for the propensity of treatment. However, the presence of dependent latent effects among the outcome and exposure models threatens the inference, since the exposure assignment mechanism is informative for the imputation of missing potential outcomes. To see this, we turn to the observed data likelihood and write

L ( θ | O ) = ∏ ν = 1 m ∫ b ν p ( Y ν , Z ν | X ν , b ν , θ ) p ( b ν | θ ) d b ν = ∏ ν m ∫ b ν ∏ i = 1 N ν p ( Y ν i | Z ν i = z ν i , Σ Z N ν i = Σ z N ν i , X ν i y , b ν y , θ ) p ( Z ν i | X ν i z , b ν z , θ ) p ( b ν | θ ) d b ν = ∏ ν m ∫ b ν ∏ i = 1 N ν p ( y ν i ( z ν i , Σ z N ν i ) | X ν i y , b ν y , θ Y ) p ( Z ν i | X ν i z , b ν z , θ Z ) p ( b ν | Σ ) d b ν

where the third equality follows from counterfactual consistency, latent exchangeability, and positivity. This means that having access to the latent confounders b would render the treatment assignment ignorable. However, since b is unobserved and yet a component of the treatment assignment mechanism through Σ , the treatment assignment is informative for the imputation of missing potential outcomes and needs to be accounted for in the Bayesian procedure.^11,13,35,36 Hence, in the presence of unobserved contextual confounding, the propensity score model is not ignorable in the inference for μ z α , μ α and functions of these quantities.

For node i in subgraph ν , y ~ ν i ( z ν i , Σ z N ν i ) acts like a missing variable that can be imputed after setting exposure to some value ( z ν i , Σ z N ν i ) and sampling from the marginal sampling distribution for the replicated potential outcomes given parameters, assigned treatments and observed outcomes

p ( y ~ ν i ( z ν i , Σ z N ν i ) | θ , O ) = ∫ b ν p ( y ~ ν i | Z ν i = z ν i , Σ Z N ν i = Σ z N ν i , b ν , θ , O ) p ( b ν | θ , O ) d b

where the equality holds under latent ignorability because of the inclusion of the propensity score model in the posterior inference for θ .

5.5. Prior specification

We adopt a combination of vague and weakly informative priors for model parameters. These parameters are the coefficients of the outcome model in (1), the coefficients of the exposure model in (3), the residual variance and the autocorrelation parameters of the SAR covariance matrix in (2), and the parameters of the random effect covariance matrix in (4).

For the coefficients of the outcome and exposure models, β and γ , we adopt vague independent N ( 0 , σ prior 2 ) prior distributions with σ prior 2 = 100 . We specify τ ∼ Uniform [ − 1 , 1 ] and σ ε ∼ Half-Cauchy ( α prior ) with the recommended hyperparameter of α prior = 25 , which results in a reasonably flat prior distribution.^37,38

For the covariance matrix of the random effects, we propose using the decomposition Σ = S R S , where S is a ( q 1 + q 2 ) × ( q 1 + q 2 ) diagonal matrix with the standard deviations of the random effects as its diagonal elements and R is a correlation matrix, and assigning the prior R ∼ LKJ ( η prior ) with η prior = 1 , where LKJ refers to the Lewandowski–Kurowicka–Joe distribution. ³⁹ This choice of prior assigns equal density to all correlation matrices of order q 1 + q 2 . Finally, we assign independent Half-Cauchy ( α prior ) priors to the diagonal elements of S . In the simulations and application that follow, we consider the joint random intercept model, that is, b ν = ( b ν y , b ν z ) , such that

S = ( σ b y 0 0 σ b z ) , R = ( 1 ρ ρ 1 )

In that case, the LKJ prior reduces to assigning a Uniform [ − 1 , 1 ] prior to ρ .

The resultant posterior distributions are approximated through computational approaches. We employ Markov chain Monte Carlo methods.

5.6. Bayesian standardisation

In this section, we describe a procedure to recover the posterior distributions of μ z α , μ α , and relevant contrasts from the posterior distribution of θ . In a point treatment setting, the Bayesian g-formula or standardisation typically entails marginalising the individual average potential outcome over the distribution of covariates, which requires obtaining the posterior distributions of θ and θ X . However, recall that in this work, we do not model the multi-dimensional covariates X and instead rely upon mixed average causal effects, which are reasonable approximations to population average causal effects. ²³

Samples from the posterior distributions of μ z α , μ α , and relevant contrasts are obtained by marginalising the posterior distribution of μ z α ( x ) ≡ μ z α ( x ; θ ) over the empirical distribution of X and the subgraph distribution, where the notation highlights the dependence on the outcome model parameter. The standardisation procedure essentially consists of computing the following integral:

μ ~ z α = ∫ ν ∫ X μ z α ν i ( x ; θ Y ) F ^ X ν ( x ) H ^ ( ν ) d x d ν = m − 1 ∑ ν = 1 m N ν − 1 ∑ i = 1 N ν μ ~ z α ν i ( x ν i ; θ )

where μ ~ z α ν i ( x ν i ; θ ) = ∑ Σ z N ν i = 0 d ν i y ~ ν i ( z , Σ z N ν i ) π ( Σ z N ν i ; α ) is the imputed individual average potential outcome and y ~ ν i ( z , Σ z N ν i ) ≡ y ~ ν i ( z , Σ z N ν i , x ν i ) is a random draw from the marginal sampling distribution for the replicated potential outcomes given individual treatment z , neighbourhood treatment Σ z N ν i , and covariate vector x ν i . We use Monte Carlo integration to obtain samples from the marginal distribution for the potential outcomes with respect to the random effects. The Bayesian standardisation procedure is summarised as follows:

Set an exposure value z and an allocation strategy α .

For each posterior draw θ ( s ) , s = 1 , … , S , perform the following steps:

For each subgraph ν ∈ { 1 , … m } , sample a random effect vector b ~ ν b ( s ) ∼ N ( 0 , Σ ( s ) ) .

For each subgraph ν ∈ { 1 , … m } , each node i ∈ { 1 , … , N ν } , and each neighbourhood exposure value Σ z N ν i ∈ { 0 , … , d ν i } , draw the corresponding missing potential outcome from the conditional posterior predictive distribution

y ~ ν i b ( s ) ( z , Σ z N ν i ) = β Z ( s ) z + h 1 ( Σ z N ν i ; β Z N ( s ) ) + h 2 ( z Σ z N ν i ; β Z × Z N ( s ) ) + h 3 ( x ν i ; β X ( s ) ) + W ν i ⊤ b ~ ν y , b ( s ) + δ ν i ( s )

Repeat steps (a) and (b) B times for each exposure value ( z , Σ z N ν i ) to obtain a draw from the sampling distribution of the missing potential outcomes conditional on θ (marginalised with respect to b ) for all nodes and all subgraphs using

y ~ ν i ( s ) ( z , Σ z N ν i ) = 1 B ∑ b = 1 B y ~ ν i b ( s ) ( z , Σ z N ν i )

Using the imputed potential outcomes in Step 2, for each posterior draw s = 1 , … , S , subgraph ν ∈ { 1 , … , m } , and node i ∈ { 1 , … , N ν } , obtain the posterior individual average potential outcome as

μ ~ z α ν i ( s ) ( x ν i ) = ∑ Σ z N ν i = 0 d ν i y ~ ν i ( s ) ( z , Σ z N ν i ) π ( Σ z N ν i ; α )

Finally, obtain the posterior population average potential outcome for treatment z and coverage α as μ ~ z α ( s ) = m − 1 ∑ ν = 1 m N ν − 1 ∑ i = 1 N ν μ ~ z α ν i ( s ) ( x ν i ) , as well as the posterior population marginal average potential outcome, μ ~ α ( s ) = α μ ~ 1 α ( s ) + ( 1 − α ) μ ~ 0 α ( s ) .

Approximations of the posterior distributions for relevant causal contrasts can be obtained by applying relevant functions ( DE ( α ) , IE ( α , α ′ ) , TE ( α , α ′ ) , OE ( α , α ′ ) ) to each posterior draw μ ~ z α ( s ) and μ ~ α ( s ) , s = 1 , … , S .

6. Simulations

We present simulation results under the joint random intercept model, where b ν = ( b ν y , b ν z ) to assess frequentist properties of the proposed Bayesian estimator based on the JMM in finite samples compared to estimators based on LMM and FEM, respectively. The LMM approach incorporates a cluster-level random intercept in the outcome model, but ignores the treatment assignment mechanism completely in the sense that the propensity score component is dropped from the likelihood altogether, while the FEM approach also drops treats the outcome data as unclustered and thus does not account for within-cluster dependence. Additionally, we compare the joint modelling approach to an augmented linear mixed-effects model approach that includes the joint propensity score as an additional covariate (LMM-augmented (LMM-aug)), which is described in more detail in Supplemental Appendix B. The default simulation is conducted according to the following steps with m = 50 disjoint subgraphs, such that A is block-diagonal.

Generation of the network: For ν ∈ { 1 , … , 30 } , we sampled N ν from the distribution Poisson ( 35 ) , and for ν ∈ { 31 , … , 50 } , we sampled N ν from the distribution Poisson ( 12 ) . For each subgraph, we generated the random vector H ν of size N ν , whose elements are from a Bernoulli ( 0.5 ) distribution. The m subgraphs were generated according to an exponential-family random graph model in which matching H values increase the probability of a tie between two nodes, inducing network homophily. The final network G was taken as the union of the m components.

For ν = 1 , … , m , we sampled ( b ν y , b ν z ) from N 2 ( 0 , Σ ) , with

Σ = ( σ b y 2 σ b y σ b z ρ σ b y σ b z ρ σ b z 2 )

and sampled the SAR errors δ ν from N N ν ( 0 , σ ε 2 U ν U ν ⊤ ) .

For ν = 1 , … , m , i = 1 , … , N ν , two baseline covariates were generated as X 1 ν i ∼ N ( 0 , 1 ) and X 2 ν i ∼ Bernoulli ( 0.5 ) . Define the proportion of treated nodes in a neighbourhood as ψ ( z N ν i ) = Σ z N ν i / d ν i . Potential outcomes for all z ν i ∈ { 0 , 1 } , Σ z N ν i ∈ { 0 , 1 , … , d ν i } were given by

y ν i ( z ν i , z N ν i ) = 2 + 2 z ν i + ψ ( z N ν i ) + z ν i ψ ( z N ν i ) − 1.5 | X 1 i | + 2 X 2 i − 3 | X 1 i | X 2 i + b ν y + δ ν i

where δ ν i is the i -th entry of the SAR errors vector δ ν .

The treatment was generated as Z ν i = Bernoulli ( p ν i ) , where

p ν i = logit − 1 [ 0.1 + 0.2 | X 1 i | + 0.2 X 2 i | X 1 i | − H ν i + b ν z ]

Under this model, nodes that share ties tended to have similar treatment assignments.

In the standard simulation scenario, simulations were carried out 500 times by repeating steps 2–4 with Var ( b ν y ) = Var ( b ν z ) = 1 and a correlation coefficient ρ = 0 , 0.2 , 0.5 , 0.8 between the outcome and exposure model random intercepts. For simplicity, we first set conditionally independent random errors for observations within each component with τ = 0 and σ ε 2 = 1 , such that the covariance matrix of the random error vector within a component was equal to σ ε 2 U ν U ν ⊤ = I N ν . The true parameters were obtained by averaging the potential outcomes that we generated in Step 2 over the simulated datasets. For each simulated dataset and each model (JMM, LMM, LMM-aug, and FEM), we fit four chains with 10,000 iterations used as a burn-in period, and thinning by 50 iterations. The posterior distributions of μ z α and μ α for z = 0 , 1 and α ∈ { 0.1 , 0.2 , … , 0.9 } as well as relevant causal contrasts were obtained from the JMM, LMM, LMM-aug, and FEM parameter posterior distributions, respectively, using the procedure detailed in Subsection 5.6 with B = 50 Monte Carlo samples. For all estimands of interest, the Bayesian estimator was defined as the posterior mean, and we calculated the posterior standard deviation and a 95% equi-tailed credible interval as measures of uncertainty.

Table 1 shows the relative bias (RB, %), the average of posterior standard deviations (ASDs), the empirical standard deviation (ESD) of Bayesian estimators, and the coverage of 95% credible intervals (empirical coverage probability (ECP, %) associated with JMM, LMM, LMM-aug, and FEM for the estimation of μ 0 , 0.7 , μ 1 , 0.7 , and DE ( 0.7 ) = μ 1 , 0.7 − μ 0 , 0.7 . Our discussion of results focuses on the direct effect DE ( 0.7 ) . In general, JMM provides valid inference in comparison to LMM, LMM-aug, or FEM, and as expected, the benefit associated with accounting for contextual confounding becomes clearer as the correlation ρ increases. The relative bias associated with LMM or LMM-aug, albeit never very large, increases with ρ , while the relative bias associated with JMM remains small. In general, the posterior standard deviations accurately estimate the variability of the regression estimators when either JMM, LMM, or LMM-aug is used, whereas the actual variability of the causal estimators is drastically underestimated when the multilevel structure of the data is not accounted for in FEM.

Figure S1, which can be found in Supplemental Appendix C, displays the average of the 500 JMM, LMM, LMM-aug, and FEM estimators of the direct effect DE ( α ) as well as the corresponding true dose–response curve, where each panel corresponds to one of the ρ values considered. Figure S2 in Supplemental Appendix C displays the average of the 500 JMM, LMM, LMM-aug, and FEM estimators of the indirect effect IE ( α , α ′ ) as well as the corresponding true surface. LMM, LMM-aug, and FEM estimators of the direct effect sustain a positive bias as soon as there is the presence of unmeasured contextual confounding ( ρ > 0 ), and this bias increases as the contextual confounding intensifies. The same conclusions pertain to the indirect effect, which takes the form of a surface as a function of α and α ′ .

In addition to the simulation scenario above, we generated datasets under a variety of conditions with a fixed ρ = 0.2 . We considered the following seven scenarios for which results can be found in Table 1:

We assumed a SAR structure for the random errors with τ = 0.1 and modelled them as such across the different methods (JMM, LMM, LMM-aug, and FEM). Comparing this scenario to the standard scenario, JMM, LMM, and LMM-aug’s performance remain largely similar.

In an additional scenario, we generated data under no treatment effect with the following outcome model

y ν i ( z ν i , z N ν i ) = 2 − 1.5 | X 1 i | + 2 X 2 i − 3 | X 1 i | X 2 i + b ν y + δ ν i

However, we included terms for the treatment, the neighbourhood treatment, and the interaction at the modelling stage. The average Bayesian estimate of DE ( 0.7 ) is closest to 0 with JMM, followed by LMM, LMM-aug, and FEM.

To mimic certain aspects of the Add Health study, we also generated the potential outcomes under a truncated normal distribution with a conditional mean

2 + 2 z ν i + p ( z N ν i ) + z ν i p ( z N ν i ) − 1.5 | X 1 i | + 2 X 2 i − 3 | X 1 i | X 2 i + b ν y

standard deviation of 1, a lower truncation bound of − 5 , and an upper truncation bound of 9. Again, we obtain comparable results to the main simulation scenario under this simulation scheme.

To assess the sensitivity of the different methods to model misspecification, we analysed generated data under the standard simulation scenario using the following incorrect outcome model:

E [ y ν i ( z ν i , z N ν i ) ] = β 0 + β 1 z ν i + β 2 p ( z N ν i ) + β 3 z ν i p ( Z N i ) + β 4 X 1 i + β 5 X 2 i

Additionally, similar to Shardell and Ferrucci, ¹⁰ we considered an additional scenario of shared random intercepts (i.e. b ν z = b ν y = b ν ∀ ν ), where the outcome and treatment assignments were generated as

y ν i ( z ν i , z N ν i ) = 2 − 1.5 | X 1 i | + 2 X 2 i − 3 | X 1 i | X 2 i + λ b ν + δ ν i

Z ν i = Bernoulli ( logit − 1 [ 0.1 + 0.2 | X 1 i | + 0.2 X 2 i | X 1 i | − H ν i + ω b ν ] )

We apply the estimators proposed in Section 5 to estimate the causal effects of college maternal education on academic performance. We investigate this question while aiming to mitigate bias from potentially missing subgraph-level confounders. In the following, the outcome is a student’s GPA, and the exposure is defined as the binary indicator of whether their mother has completed four years of college. The students who did not report living with their biological mother, stepmother, foster mother, or adoptive mother at the time of the survey were excluded from the analytic sample. We also excluded isolates from the sample for computational reasons. There was partial actor non-response in the dataset, which occurs when some, but not all, attribute information is available for one or more nodes in the network. ⁴¹ To remedy this, a single imputation of the nodal attributes was performed using imputation by chained equations with predictive mean matching and 20 cycles. The second column of Table S5 in Supplemental Appendix D provides descriptive statistics for the nodal attributes after the imputation. The network structure of the Add Health subsample of m = 16 schools had 29 connected components with 7,931 students and 26,384 shared connections (average degree was 6.65) after the exclusions. Among the 7,931 students, 2,873 (36.2%) had a college-educated mother. Figure S3 in Supplemental Appendix D displays a histogram of the observed proportion of friends with a college education in the analytical sample, which shows a right-skewed distribution. The following baseline covariates were included in the outcome model: sex, age, student’s race (as a proxy for the race of the mother), whether the student was adopted, whether the mother was born in the US, whether the father is at home, household size, whether the mother works for pay, student’s perception of how much their mother cares for them, club participation, motivation at school, trouble at school, sense of belonging, attendance, health status, physical fitness, and screen time. The model for the exposure only included four of the baseline covariates: the student’s race (as a proxy for the race of the mother), whether the mother was born in the US, whether the father is at home, and household size. Checks for positivity were performed by comparing the distribution of the individual propensity score for the exposure model, including the four potential confounders across the two groups defined by maternal education status. Figure S4 in Supplemental Appendix D shows reasonable covariate overlap across the two exposure groups.

We considered three analyses (JMM, LMM, and FEM), all of which adjust for the same baseline covariates enumerated above and include all variables linearly in the predictors of the outcome model (and of the propensity score model in the case of JMM). To conform with the observed support of the GPA distribution, we considered a truncated Gaussian likelihood for the outcome, with the lower and upper truncation bounds set to 1 and 4, respectively. The main analyses do not impose a SAR structure on the random errors, since the adjacency matrix in the Add Health subnetwork is not block-diagonal (there are ties between students from different schools) and its inversion would be prohibitively computationally expensive considering its large dimension ( N = 7 , 931 ).

Estimates based on the JMM, LMM, and FEM approaches for DE ( α ) , IE ( α , 0.3 ) , TE ( α , 0.3 ) and OE ( α , 0.3 ) along with 95% credible intervals are shown in Figure 3 and Table S6 in Supplemental Appendix D. In general, JMM, LMM, and FEM estimators give similar conclusions, in that the directions of the dose–response curves align. The point estimates suggest that academic performance is increased not only directly by having a college-educated mother, but also indirectly by having a greater proportion of friends with a college-educated mother. We provide an interpretation for DE ( 0.2 ) . If every student in this network were to be “assigned” a college-educated mother versus a mother with a non-college education, we would expect the average GPA difference to be 0.0252 (95% CrI: [-0.0291, 0.0792]) points under a counterfactual probability of receiving treatment α = 20 %. Using the JMM method, we obtain an estimate of 0.0025 (95% CrI: [-0.0373, 0.0428]) for IE ( 0.8 , 0.3 ) . This implies that if we held the maternal education indicator to 0, the GPA of students with 80% of friends with college-educated mothers would be on average higher by 0.0025 points compared to a setting in which the maternal education indicator setting is again set to 0 with, however, 30% of friends having college-educated mothers, assuming both groups’ covariate distribution matched that of the overall population.

OPEN IN VIEWER

Figure 3.

JMM, LMM, and FEM estimates of DE ( α ) , IE ( α , 0.3 ) , TE ( α , 0.3 ) , and OE ( α , 0.3 ) effects with 95% credible intervals. JMM: joint mixed-effects model; LMM: linear mixed model; FEM: fixed effects model.

Most credible intervals shown in Supplemental Table S6 contain zero, with the exception of DE ( 0.4 ) and TE ( α , 0.3 ) for α = 0.4 , 0.6 , 0.8 when FEM is used. In general, effect estimates are closer to zero based on JMM and LMM compared to the FEM approach, which ignores the clustering in the graph. Effect estimates based on our approach (JMM) are also lower in magnitude than those based on LMM, which assumes ignorability of the exposure assignment. In the JMM approach, the posterior median for ρ was 0.58 (95% CrI: [ 0.09 , 0.84 ] ). These results suggest that there might exist unmeasured contextual confounding, and our approach partially mitigates the bias due to the unmeasured latent traits. Also, the credible intervals for IE ( α , 0.3 ) based on JMM are slightly narrower than those for LMM, illustrating that accounting for unmeasured contextual confounding can also lead to efficiency gains in the estimation of interference effects.

In this paper, we developed a Bayesian approach for causal inference with clustered network-based observational studies, estimating the direct and indirect effects of a binary exposure on a continuous outcome while accounting for potential unmeasured contextual confounding. We discussed the inherent challenges that arise in causal inference with multilevel network data, which are prevalent in education, and discussed the relevance of accounting for potential non-ignorability of the treatment assignment. We proposed to simultaneously model the outcome and treatment assignment by JMM, which reduces unmeasured contextual confounding when models are correctly specified. Our proposed methodology can also address potential network autocorrelation in the outcome. We described a Bayesian standardisation procedure that can be used to obtain the posterior distributions of causal estimands of interest from the posterior distributions of the JMM parameters. We note that JMM has a specific representation of unmeasured contextual confounding, the level of which is governed by the level of heterogeneity in the outcome and treatment assignment mechanisms and the correlation between the treatment assignment and the outcome. In the case of the joint random intercept model, this correlation operates through the correlation parameter ρ . Large random effect variances or ρ values may lead to significant bias due to unmeasured confounding. The simulation study demonstrated that JMM provides valid inference in a wide variety of settings for average potential outcomes and causal contrasts, where LMM and FEM may not. A particular advantage of our proposed approach is that it performs well with either the presence or absence of unmeasured contextual confounding.

Our proposal fits more widely into the literature of Bayesian causal inference methods, incorporating a model for the propensity score. Existing ways of combining propensity score and outcomes models include the propensity score as a covariate in the outcome model, specifying dependent priors or shared parameters between the propensity score and outcome models, or posterior predictive inference through inverse probability weighting and doubly robust estimation.^23,42,43 As we used a latent Gaussian structure to relate the outcome and propensity score models, our method falls under the second category of methods. ¹³

As previously mentioned, JMM relies on parametric assumptions regarding the joint distribution of outcomes and treatment assignments. Unbiased estimation requires correct specification of the outcome model, the treatment assignment mechanism, and the random effects distribution. Under the assumption of strong ignorability, parametric model assumptions would not necessarily be required for valid causal inference. However, in the context of a school-based design such as Add Health, unmeasured school-level confounders are highly likely, as students attending the same school may share latent characteristics affecting both their outcome and exposure. To relax the unconfoundedness assumption, the trade-off here is to make additional assumptions on model specifications. In particular, we assumed that random effects follow a normal distribution in both the outcome model and the propensity score model. As seen in simulation studies, a violation of the Gaussian assumption in the form of exponentially distributed random effects can impact the coverage probability of credible intervals for average potential outcomes more severely, while inference for causal contrasts seems to be mildly affected for the data-generating mechanism considered. Importantly, similar to LMM and LMM-aug, JMM also yields biased inference under outcome model misspecification, which aligns with theoretical expectations. As it is more restrictive and also computationally expensive, the gain provided by JMM may not be worthwhile in larger networks or when contextual confounding is weak. Another caveat specific to our analysis of spillover effects of maternal education in Add Health is that we had to use a smaller subset of the dataset for computational reasons, potentially leading to an underpowered analysis. In addition, we considered a linear functional form to relate the outcome to the predictors in our analysis of Add Health. While posterior predictive checks and the deviance information criterion supported the use of a linear model over a generalised additive model in this particular case, in general, non-linear such as spline representations, as the general outcome model formulation in (1) should be investigated, as failure to capture true non-linear relationships may lead to invalid causal inferences.

Despite the potential merits of this work, methodological challenges remain in order to improve inference for causal effects in the presence of network interference. We performed approximate checks for positivity by assessing covariate overlap for the individual exposure only. However, the context of interference requires sufficient overlap across groups defined by the different combinations of the individual and neighbourhood exposures. Designing positivity assessments for the network interference setting should be included in future research. Also, the impact of the network topology on the covariance matrix Φ in the SAR model remains to be investigated. For instance, a sparse or inhomogeneous graph may have an effect with regard to the computation of Φ , particularly on the matrix inversions involved. Although the methodology in this paper is developed under a SAR structure, the analysis of the Add Health study was done under τ = 0 corresponding to independent random errors, avoiding the SAR structure for this large-scale empirical network of N = 7 , 931 nodes. In theory, for large sparse networks, one can avoid inversion entirely by parameterising the model in terms of the precision matrix rather than the covariance matrix. This approach, in combination with an incomplete Cholesky factorisation, would allow for scalable inference without ever forming or inverting the full covariance matrix. While our implementation relied on the covariance form, a precision-based parameterisation is a natural avenue for handling larger networks in future work.

Of the different sources of unmeasured network confounding, ⁶ our approach only addresses contextual confounding. However, node-level unmeasured confounders could still introduce bias in causal effect estimates. One direction for future research could entail developing novel sensitivity analysis techniques to assess the impact of node-level unmeasured confounding on the inference for direct and spillover effects, such as a method based on the confounding method, which has the advantage of avoiding introducing a hypothetical unmeasured confounder.^44,45

Another challenge inherent to this context is homophily bias. In sociological settings, this can lead to latent variable dependence due to latent traits that both drive the formation of network ties and predict an individual’s treatment or outcome. In all scenarios we considered here, we assumed that all variables driving homophily were measured and conditioned upon in analyses, which is unlikely to be the case in practice. Imposing a SAR structure on the random errors of the outcome model might account for some of the latent dependence, but more work has to be done to better understand how homophily can affect inference in these settings. Latent homophily poses a threat to the identification of causal effects in various ways. Shalizi and Thomas ⁴⁶ showed that unless all of the nodal attributes driving both social-tie formation and the behaviour of interest are observed, then peer effects are generally unidentified. Assuming that the social network is generated by a latent community, McFowland III and Shalizi ⁴⁷ proposed asymptotically unbiased estimators of peer influence under homophily by controlling the latent location of each node. Another avenue for future work would be to use their methodology and incorporate estimates of latent positions in the outcome and exposure models to control for potential latent homophily when inferring for interference effects.

Another limitation of this work is that it is assumed that the adjacency matrix A is fully observed and without error. A common strategy consists of treating missing edges as absent or dropping vertices with missing edge information, which can lead to dubious results.^48–50 Li et al. ⁵¹ have studied the impact of noisy networks for experiments on networks, but the problem of missing data remains largely understudied in the network interference literature. To our knowledge, the only example of a use of network imputation to assess causal effects with interference is in the context of experiments on networks. ⁵² This could be another important avenue for future work in the context of network-based observational studies.

In general, cross-sectional social network data are limited, since it remains a challenge to distinguish whether associations in peers’ outcomes are due to interference or homophily. ⁶ Future developments should entail extending methods presented in this paper to discrete-time dynamic networks and considering the co-evolution of network ties, A i j t , and vertex attributes, ( Y ν i t , Z ν i t , X ν i t ) , through time t for more nuanced causal inferences. A more straightforward path would be to consider discrete-time processes, as discrete-time models for time-varying effects, such as structural nested models or marginal structural models estimated via the g-formula, g-estimation, or related methods, do not generalise easily to continuous time, due to the entanglement of uncountably infinite variables. Stochastic actor-oriented models are appealing in that they allow for joint modelling of social influence and social ties, and can hence provide adequate control for network endogeneity in longitudinal social network data.^6,53 Despite their potential for interference analyses, these models have yet to be formulated in a counterfactual framework. ⁶ Special considerations would be needed to account for potential time-varying confounding by X ν i t as well as by A i j t , since these variables could also be on causal paths for subsequent time points. This represents a promising yet challenging avenue for future work. Lastly, future extensions of the JMM could consider categorical and count outcomes as well as survival outcomes.

Supplemental Material