2.1 Modeling interventions with LSMs

Consider an intervention in which entire networks are assigned to treatment or condition. Because we will eventually model the effect of the intervention as a function of $\sigma_Z^2$ in Equation (1), we first demonstrate how $\sigma_{Z}^2$ impacts network structure.

Recall that in a LSM, each network is represented by relative positions of its nodes in a low-dimensional space and the probability of a tie is a function of the distance between the nodes. Thus, $\sigma_{Z}^2$ controls the spread or scale of the distance, and networks with smaller values of $\sigma_{Z}^2$ will have higher density than networks with larger values of $\sigma_{Z}^2$ , because smaller values of $\sigma_{Z}^2$ map to smaller pairwise distances and higher tie probabilities. In addition, we also expect to see an increase in the number of isolated nodes as $\sigma_{Z}^2$ increases.

However, unlike random graph models, LSMs are also known to preserve reciprocity and transitivity (Hoff et al., Reference Hoff, Raftery and Handcock2002). This is because, as long as the position of a node relative to the neighboring nodes is preserved, reciprocity and transitivity will also be preserved. Therefore, with the exception of extreme differences in $\sigma_{Z}^2$ , reciprocity and transitivity should be similar across values of $\sigma_{Z}^2$ .

For illustration, we simulated 100 networks for values of $\sigma_{Z}^2$ and plotted the median network statistic value (orange circle) along with 0.025 and 0.975 quantiles, indicated by the line segment. Figure 1 illustrates how the variance of the latent space positions in a LSM impacts density, variance in degree, the number of isolated nodes, and mean geodesic distance (among connected nodes). The plot showing geodesic distance is misleading. The distances are not actually decreasing with increasing $\sigma_{Z}^2$ ; rather, the number of isolated nodes who are not connected to other nodes is increasing, and these distances are not included.

Figure 1. Network statistics from networks generated from latent space models where $\sigma_Z^2$ ranges from 2 to 10 to illustrate the impact of $\sigma_Z^2$ on density and the number of isolated nodes with comparatively little impact on reciprocity and transitivity.

Aside for extreme differences in $\sigma_{Z}^2$ , there is little impact on the levels of reciprocity and transitivity although we note that reciprocity and transitivity naturally increase as density increases. Here, we use $\sigma_{Z}^2$ values 2 to 10 to illustrate a wide range, but note that in Section 4.1, we consider networks generated from LSMs with expected values of $\sigma_{Z}^2$ of 2.7 in the control condition.

Thus, an intervention that impacts $\sigma_{Z}^2$ would change the overall density and the number of isolated individuals in the network, while preserving some standard network structures such as transitivity and reciprocity. For example, in an intervention focused on improving school climate, we may expect greater unity and consensus among its teachers, but we would not expect clique structure or reciprocity to change. Similarly, teachers who were isolated previously may be now be connected to someone in the network as a result of a decrease in $\sigma_{Z}^2$ .

3.1 Estimation

The HLSMMs in this manuscript were fit using a Markov chain Monte Carlo (MCMC) algorithm (Gelman et al., Reference Gelman, Carlin, Stern and Rubin2013), and we used a combination of Metropolis-Hastings and Metropolis updates within Gibbs sampling. The Metropolis/Gibbs algorithm is standard and allows any reader familiar with MCMC to understand the estimation process as well as augment the code to expand or adapt the proposed model for their purposes. However, the HLSMM can be fit using a variety of MCMC algorithms as well as posterior approximation methods such as variational inference.

Moreover, the MCMC algorithm combines the network portion of the model and the mediation model and updates all parameters simultaneously. That is, the latent positions given the observed networks are updated iteratively alongside the regression coefficients of the mediation portion of the model. Fitting the HLSMM in this way accounts for the estimation error in the latent positions when estimating the regression coefficients used to calculate the indirect effect.

Given the normal distributions specified in the model for $\log(\sigma_{Z})$ and $\log(\sigma_{Y})$ , we used Gibbs updates for $\gamma_0$ , $\gamma_1$ , $\theta_0$ , $\theta_1$ , and $\nu$ , using normal distributions as conjugate priors. Similarly, for variance hyper-parameters $\tau_{Z}$ and $\tau_Y$ , we employed conjugate inverse-Gamma distributions so that they were also updated via Gibbs. The latent space positions as well as the variance of the latent space positions (reparametrized as $\log(\sigma_Z)$ ) were each updated using a random walk Metropolis-Hastings update. We used conjugate priors when possible, but other prior distributions can be used. For example, half-Cauchy distributions could replace inverse-Gamma distributions for the variance updates.

Note that LSMs have two identifiability issues. The first is that there are an infinite number of latent space position configurations possible due to translation, reflection, and rotation. A procrustes transformation as used in Hoff et al. (Reference Hoff, Raftery and Handcock2002) can be applied to obtain estimates of the latent space positions with a stable configuration. Note that this transformation can be applied post-sampling as convergence can be obtained without it. The second identifiability issue is discussed less often but has larger impacts on the HLSMM parameters. There is an identifiability between $\beta_0$ , which can be thought of as an intercept and related to the overall probability of a tie, and the variance of the latent space positions $\sigma^2_{Z_k}$ , which also impacts tie probability through the dispersion of nodes in the latent space. We recommend fixing $\beta_0$ when estimating $\sigma^2_{Z_k}$ is of particular importance.

4.1 Simulation study

The purpose of this small simulation study is to examine the influence of various parameters on both parameter recovery and model utility. We consider nine different combinations of parameters $\gamma_1$ and $\theta_1$ in Equation (3). Recall that $\gamma_1$ is the effect of the intervention T on the log standard deviation of the latent space positions $\log(\sigma_{Z})$ , so we explore settings where the intervention has three different effects on this measure of variability. Similarly, $\theta_1$ is the conditional effect of $\log(\sigma_{Z})$ on $\log(\sigma_{Y})$ given T, and we consider three different values for this parameter as well.

We selected values of $\gamma_1$ that resulted in differences in $\log(\sigma_Z)$ that produced visually distinct but still realistic networks and values of $\theta_1$ that yielded analogous values for $\log(\sigma_Y)$ . Note that negative values of $\gamma_1$ result in smaller variance of latent space positions. Our positive value of $\theta_1$ maps to a positive relationship between latent space variability $\log(\sigma_Z)$ and variability in the outcome $\log(\sigma_Y)$ .

Table 1 summarizes the nine cells of our simulation, and each cell represents one combination of values for $\gamma_1$ and $\theta_1$ . We also include sample variances (across replications) to illustrate the indirect and direct effects on the outcome; note that $\sigma^2_{Z}$ and $\sigma^2_{Y}$ are given instead of $ \log(\sigma_{Y})$ and $\log(\sigma_{Y})$ for ease of interpretation. The full data generating model used is given as:

(4) \begin{align}\log(\sigma_{Z_k}) &\sim N(0.5+\gamma_1T_k,~0.09), k=1,\dots,30\nonumber\\[5pt]Z_{ik} &\sim MVN_2(0, \sigma_{Z_k}^2{I}_2), i=1,\dots, 20\nonumber\\[5pt]p_{ij} &\sim \frac{\exp(-\vert Z_{ik}-Z_{jk}\vert)}{1+\exp(- \vert Z_{ik}-Z_{jk}\vert )}, i,j=1,\dots, 20\end{align}

\begin{align}{A}_{ijk} &\sim \mbox{Bernoulli}(p_{ij}), i,j=1,\dots, 20\nonumber\\[5pt]\log(\sigma_{Y_k}) &\sim N(-0.45+\theta_1 \log(\sigma_{Z_k})-0.1T_k,~0.09)\nonumber\\[5pt]Y_{ik}&\sim N (0, \sigma^2_{Y_k})~\nonumber\end{align}

where k indexes the 30 networks each with 20 nodes. The independent variable $T_k$ indicates whether network k is in the treated group such that 15 networks are treated and 15 networks are not.

Table 1. A summary of parameter values across the nine cells of the simulation study along with approximate sample variances

In Equation (4), the expected $\log(\sigma_Z)$ for control networks and treated networks is defined as $\gamma_0$ and $\gamma_0+ \gamma_1$ , respectively. Thus, in selecting values for each, our goal was to generate networks that reflected real-world data. The expected $\sigma^2_Z$ under the control condition is approximately 2.72 and under the three treated conditions ( $\gamma_1 = -0.2, -0.4, -0.6$ ) is 1.82, 1.22, and 0.82 respectively.

Figure 3 shows 30 networks generated from Equation (4) where $\gamma_0=0.5$ and $\gamma_1=-0.4$ . The control networks are shown on the top three rows of plots and the treated networks on the bottom three rows. Overall, the control networks tend to have more isolated nodes and are less dense than the treated networks. Thus, $\gamma_1=-0.2$ would result in less noticeable differences between the conditions, and $\gamma_1=-0.6$ would result in more noticeable differences between conditions.

Figure 3. Networks simulated from a HLSM for Mediation where T indicated treatment and the treated networks have smaller $\sigma_{Z_k}$ (see Equation (4)).

We selected values for $\theta_1$ in a similar manner aiming for $\sigma_Y$ to be distinct enough to detect differences. Again, we selected error variances ( $\tau=0.09$ ) to observe adequate separation of distributions of $\log(\sigma_{Y_k})$ to detect effects. Finally, we also chose to include a nonzero value for $\nu=-0.1$ to indicate a partial mediation. Figure 4 shows a plot of $\log(\sigma_{Y_k})$ against $\log(\sigma_{Z_k})$ for the 30 networks generated from Equation (4) where $\gamma_0=0.5$ , $\gamma_1=-0.4$ , $\theta_0=-0.45$ , $\theta_1=0.9$ , and $\nu=0.1$ . There is evidence of a positive effect of $\log(\sigma_{Z})$ on $\log(\sigma_{Y})$ .

Figure 4. The scatterplot of $\log(\sigma_{Y})$ against $\log(\sigma_{Z})$ shows a positive effect of the mediator on the outcome.

Within each simulation condition given in Table 1, we simulated data from a HLSMM generating model, Equation (4), and fit a HLSMM, Equation (3), using the following priors:

(5) \begin{align}\begin{split}&\gamma_0, \theta_0 \sim N(0,100)\\[5pt]&\gamma_1, \theta_1, \nu \sim N(0,9)\\[5pt]&\tau_Z, \tau_Y \sim IG(5, 0.5)\end{split}\end{align}

We selected weakly informative priors for $\gamma_0$ and $\theta_0$ and chose more informative priors for $\gamma_1$ , $\theta_1$ , and $\nu$ and error variances $\tau_Z$ and $\tau_W$ . (See the Appendix for a prior sensitivity analysis.) For each model fit, we ran three chains of length 10,000; we removed the first 1,000 draws and thinned each chain by 35 to obtain a final posterior sample of 774; this length and thinning provided adequate convergence and autocorrelation across all parameters. This process was repeated 100 times for each cell of the simulation (Table 1).

To evaluate model performance under each setting, we report the coverage rate, defined as the proportion of replications in which the 95% equal-tailed credible interval covered the data generating parameter. In addition, we also report the posterior probabilities that the mediation effect $\gamma_1\theta_1<0$ . Note that the mediation effect is negative for all cells of the simulation.

Figure 5 summarizes the posterior distributions of the mediation effect (indirect effect) $\gamma_1\theta_1$ of the 100 replications across all 9 conditions. The circles represent the maximum a posteriori (MAP) estimates, and the bars indicate the 95% highest posterior density (HPD) credible intervals. The coverage rates for the mediation effect are between 0.91 and 0.98, indicating that the indirect effect coverage is as expected. In addition, parameter coverage for all parameters estimated in the model range from 0.91 to 1.00. Furthermore, we note that the posterior variance for $\gamma_1\theta_1$ increases as the effect sizes of $\gamma_1$ and $\theta_1$ increase.

Figure 5. Simulation summary of mediation effect corresponding to settings 1–9 in Table 1. The lower and upper bound of each error bar indicates the 95% highest posterior density credible interval for the posterior distribution of the average causal mediation effect and the circle represents MAP of the distribution. The red horizontal line represents the true value of the parameters used to generate the data.

We also examine the posterior probability of the indirect effect being less than 0. We find that as the true indirect effect becomes more negative (a stronger effect), the posterior probability increases. It appears that changing the effect sizes of $\gamma_1$ has a larger impact than $\theta_1$ on these posterior probabilities. The exact probabilities are given in Table 2.

Table 2. Mean posterior probability that $\gamma_1\theta_1<0$ across all nine simulation conditions

4.2 Application to teacher networks

We now present an application of the HLSMM to real-world data. These data come from a longitudinal study of workplace and social interactions among school staff that are part of the Distributed Leadership Studies at Northwestern University (distributedleadership.org). These data represent language arts advice-seeking relationships among school staff in 14 elementary schools in 2010 and 2015 in a suburban school district in the midwestern part of the US. There were major policy changes during this time including a new mathematics curriculum, the use of mathematics instructional coaches for teachers, and high stakes testing in mathematics.

While advice-seeking ties around mathematics have been studied in detail (Spillane et al., Reference Spillane, Shirrell and Sweet2017; Spillane et al., Reference Spillane, Hopkins and Sweet2015; Spillane et al., Reference Spillane, Hopkins and Sweet2018), language arts advice-seeking ties have not. We aim to examine the effects of the policy changes in mathematics on the language arts advice-seeking networks, and whether these networks then mediated consensus about teachers’ influence on school policy.

School staff in these 14 schools were asked “During this school year, to whom have you turned for advice and/or information about curriculum, teaching, and student learning?” For each person nominated, respondents then classified the connection as being about language arts, math, or science. The language arts networks for the years 2010 and 2015 are shown in Figures 6 and 7, respectively. Note that school staff includes teachers and administrators as well as instructional support teachers (math coaches) but does not include all staff employed by the school district such as nurses, administrative assistants, custodians, or para-educators.

Figure 6. Language arts advice-seeking networks among school staff in 14 elementary schools in 2010; note that response rates in one school were particularly low.

Figure 7. Language arts advice-seeking networks among school staff in 14 elementary schools in 2015.

Response rates varied by school but on average were approximately 85%. Respondents who did not submit a survey were excluded from the network unless they were nominated by at least two other staff members, and these individuals were assumed to have no outgoing ties. Further, staff members who took the survey but did not nominate anyone remained in the networks as having no outgoing ties.

As part of this study, school staff were asked a number of items about mathematics and mathematics pedagogy, but there were a few survey items that were not subject-specific. One is a series of items asking respondents “How much influence do teachers have over school policy in each of the areas below?” and the areas listed were as follows: hiring professional staff; planning how discretionary funds should be used; determining which books and instructional materials are used in classrooms; establishing the curriculum and instruction program; determining the content of in-service programs; setting standards for student behavior; and determining goals for improving the school.

To demonstrate how the HLSMM could be used in practice, we examine the extent to which language arts advice-seeking networks mediated the effect of the district changes on consensus about teacher’s influence on school policy. The belief is that structural differences in the language arts networks in 2010 compared to 2015 may impact differences in school staff consensus about teacher influence over policy (in 2010 compared to 2015).

The HLSMM fit to these data is given as:

(6) \begin{align}A_{ijk} &\sim \mbox{Bernoulli }(p_{ij}), i,j=1,\dots, n_k\nonumber\\[4pt]\mbox{logit}\, (p_{ij}) &\sim \beta_0- \vert Z_{ik}-Z_{jk}\vert \nonumber\\[4pt]Z_{ik} &\sim MVN_2({0}, \sigma_{Z_k}^2{I}_2), i=1,\dots, n_k\nonumber\\[4pt]\log(\sigma_{Z_k}) &\sim N(\gamma_0+\gamma_1 T_k,\tau_Z)\end{align}

\begin{align}\log(\sigma_{Y_k}) &\sim N(\theta_0+\theta_1 \log(\sigma_{Z_k})+\nu T_k,\tau_Y)\nonumber\\[4pt]\gamma_0, \theta_0 &\sim N(0,100) \nonumber\\[4pt]\gamma_1, \theta_1, &\nu \sim N(0,9) \nonumber\\[4pt]\tau_Z, \tau_Y &\sim \text{IG}(5, 0.5) \nonumber\end{align}

where $A_{ijk}=1$ if teacher i seeks language arts advice from teacher j in school k and 0 otherwise, $T_k=1$ if the network is from 2015 and 0 otherwise, and $\sigma_{Y_{k}}$ is the standard deviation in teacher influence scores for school k.

Because our model hinges on accurate estimation of $\sigma_{Z_k}$ , we chose to fix $\beta_0$ . To fit the HLSMM, we ran a single chain of length 200,000. We broke the chain into three pieces to assess convergence; the potential scale reduction factor for the mediation effect was 1 with an upper confidence bound of 1.01, indicating adequate mixing. After tossing the first 1000 steps and thinning by every 175th step, our posterior sample size was 1138.

The posterior distributions for $\gamma_1$ , $\theta_1$ , $\nu$ , and $\gamma_1\theta_1$ are shown in Figure 8, where the HPD credible intervals are indicated by the shaded region. There is evidence of an intervention effect $\gamma_1$ on the $\log(\sigma_Z)$ , and the posterior probability that $\gamma_1>0$ is 0.94. The evidence of a mediated or indirect effect of $\log(\sigma_Z)$ is weak; the posterior probability that $\gamma_1\theta_1<0$ is 0.71, which is not statistically significant from a frequentist perspective that offers little evidence of a mediated effect. Further, the posterior probability that $\nu<0$ is 0.93 which suggests some evidence of a direct effect.

Figure 8. Density plots of the posterior distributions from HLSMM fit to teacher advice-seeking data. The distribution of the mediation effect $\gamma_1 \theta_1$ suggests very little evidence that the network structure mediates changes in teacher consensus.

Thus, we conclude that there is evidence that the networks in 2015 compared to those in 2010 have a larger latent space variance in a LSM. There is some evidence of a direct effect of policy changes on staff consensus about their influence, but little evidence that the spread of the latent space positions mediated the effect of policy changes on teacher consensus.

This example serves as an illustration of our model as opposed to a formal analysis of these networks. Our evidence was quite weak, but even if our finding for a mediated effect could be considered statistically significant, there were too many limitations for these results to be used in practice. We are comparing networks in 2010 with networks in 2015, and there may be other causes for differences in network structure over the years than just the district-led policy changes. Similarly, there may be other variables that impact differences in teacher consensus.

As an illustration, however, we can continue with this real-world data analysis to examine goodness of fit. Note that we did fix $\beta_0$ as there is potentially an identifiability issue between $\beta_0$ and $\sigma_{Z_k}$ . We selected several possible values for $\beta_0 = (-3,-2,-1, 0, 1, 2,3)$ . If we fixed the intercept at a value that was too small, the overall probability of tie would be not be well recovered because the model is defined by subtracting distance between the latent positions from the intercept. Subtracting a number from a small number would yield results in fitted model that predicts extremely few ties. Similarly, if we selected an intercept that was too large, the latent positions would need to be quite far and it is possible that our model would produce networks dissimilar to our observed data. We hypothesized that $\beta_0=1$ would be ideal and decided to rely on goodness-of-fit checks to confirm.

We use posterior predictive checks and generate network data from the parameters specified at each step in our MCMC chain. We examine goodness of fit in two stages: we first investigate the network portion of the model and then we examine the regression portion of the model. For the network model fit, we use the latent space positions estimated at each step of the chain to generate 28 networks and evaluate our model fit using density, reciprocity, and transitivity in each of the 28 networks.

Figure 9 shows the distributions of density, reciprocity, and transitivity for each network generated from parameters at each step in the MCMC chain. The observed network statistic for each network are shown with a vertical line. In general, density, reciprocity, and transitivity are all well recovered, indicating that the LSM was able to capture these aspects of the observed network.

Figure 9. Posterior predictive checks for the latent space portion of the HLSMM show that network density, reciprocity, and transitivity are well recovered.

Note that posterior predictive checks using values of $\beta_0 =-1, -2$ resulted in models that generated networks with lower densities, reciprocity, and transitivity values than in our observed network; similarly values of $\beta_0=3$ produced networks that tended to generate networks higher in reciprocity and transitivity than in our observed data. Using $\beta_0=0$ or $\beta_0=2$ did not produce extremely different networks, although reciprocity recovery was slightly lower for $\beta_0=0$ .

For the regression portion of the model fit, we took the parameters at each step of the MCMC chain to generate outcome data $\log \sigma_{Y_k}$ . We show the distribution of our predicted outcome data with the observed data in Figure 10. The posterior predictive checks suggest the outcome data were well recovered.

Figure 10. Posterior predictive checks for the regression portion of the HLSMM show that the predicted outcome $\log(\sigma_Y)$ in each network is well recovered.

Thus, we find the HLSMM fitted to the language arts advice-seeking network data to be adequate because the HLSMM generates data quite similar to our observed data. Despite the lack of evidence of mediation, we believe this example illustrates the feasibility of our model and demonstrates how posterior predictive checks can be used to assess model fit.