A Bayesian multifactor spatio-temporal model for estimating time-varying network interdependence

Abstract

This paper proposes a Bayesian multilevel spatio-temporal model with a time-varying spatial autoregressive coefficient to estimate temporally heterogeneous network interdependence. To tackle the classic reflection problem, we use multiple factors to control for confounding caused by latent homophily and common exposures. We develop a Markov Chain Monte Carlo algorithm to estimate parameters and adopt Bayesian shrinkage to determine the number of factors. Tests on simulated and empirical data show that the proposed model improves identification of network interdependence and is robust to misspecification. Our method is applicable to various types of networks and provides a simpler and more flexible alternative to coevolution models.

Network influence refers to the intervening effect of a social network on actors’ choices and behavior, which leads to correlated nodal outcomes. In today's hyper-connected world, network influence is so ubiquitous and important that it is no longer a nuisance parameter but has become a quantity of primary interest in many substantive and methodological studies. This research focuses on identifying and estimating network interdependence, which is also known as endogenous network effects.¹ While exogenous network effects measure how alters’ attributes or previous outcomes affect ego's outcome through network channels, network interdependence refers to contemporaneous and reciprocal effects between nodal outcomes. When studying network interdependence, social networks are conceptualized as payoff structures and strategic settings in which actors’ interactive choices form an equilibrium system (Elkins and Simmons, 2005), and strategic interaction and optimizing behavior of networked actors are deemed as the sources of outcome interdependence (Manski, 1993; Franzese and Hays, 2008b). Network interdependence is substantively interesting because it generates social multipliers and has unique implications for policy interventions (Burke, 2016).

This paper proposes a multilevel spatio-temporal model with a multifactor error structure (henceforth, MLST-MF) to estimate time-varying network influence with longitudinal data. The spatial autoregressive regressions (SARs) are the most employed models to estimate network interdependence because SARs consider first- and higher-order connectivity and the spatial autoregressive coefficient (ρ) estimates the quantity of interest (Bramoulle et al., 2014). With fast accumulations of network data, studies on network interdependence are often longitudinal, but the existing spatio-temporal models mostly assume a constant ρ. With the dimension of time, we can no longer assume interdependence to be invariant when networks change over time and so do their implications to social interaction. The proposed model utilizes the multilevel setting to allow the spatial autoregressive coefficient to vary over time and include macro-level covariates (e.g., structural features of the network) to explain the variation of network interdependence. This specification suggests network interdependence is jointly determined by micro-level behavior and macro-level structural features.

The second important feature of the proposed model is the multifactor error structure. The latent factor approach has been widely applied in the causal inference literature as a bias-reduction technique (e.g., Xu, 2017; Samartsidis et al., 2020; Pang et al., 2021). We adopt the approach as a main strategy to improve identification of network interdependence. Identification requires to separate network influence from other sources of correlated outcomes, such as homophily and common exposures, which is known as the reflection problem (Manski, 1993) and the inverse Galton's problem (Gleditsch and Ward, 2006). The challenge can be illustrated with the directed acyclic graph approach (Pearl, 2009): consider a backdoor path, y _jt ← {z _j,  z _j,t−p} → w _ijt ← {z _i,  z _i,t−p} → y _it, where y _it and y _jt denote outcomes of ego and alter, w _ijt is a network tie from ego i to alter j at time t, z's are time-invariant or time-varying attributes of ego and alter, and p is a non-negative integer. It is easy to see that w _ijt is a collider in this path and should not be included in the estimation of the relationship between the outcomes. But by definition network effect analysis inevitably involves w _ijt, which opens the backdoor path and invites endogeneity bias (Shalizi and Thomas, 2011; Elwert and Winship, 2014). To identify network interdependence we need to re-block this path by controlling for either {z _j,  z _j,t−p} or {z _i,  z _i,t−p}. But if there are latent nodal or contextual attributes in these sets, the backdoor path is still open and the estimate of network interdependence is biased.

Because the identification problem is rooted in the inter-connection of selection and influence, an ideal solution is to endogenize both processes. However, the behavior equation in coevolution models, such as extended stochastic-actor-oriented models (SAOMs), does not have the interpretation of an equilibrium system (Han et al., 2021), and endogenous network influence is seldom the quantity of interest in such models. More importantly, coevolution models are still subject to latent confounders (Robins, 2015; Ragan et al., 2019; He and Hoff, 2019). Furthermore, they impose restrictions on network ties and are applicable primarily to networks built on binary ties (Snijders, 2001, 2005; Steglich et al., 2010). Because of these limitations and restrictions, researchers have developed and applied alternative methods with bias-correction techniques and simpler model specifications to identify network influence. For instance, Christakisa and Fowler propose a model including ego's lagged outcome, arguing that latent confounders affect current outcomes mostly via previous outcomes (Christakisa and Fowler, 2008, 2013). This argument has been widely challenged and their model also suffers from other estimation and inference problems (Robins, 2015). Another popular approach in applied studies is to use two-way fixed effects. Two-way fixed effects can only capture unit- or time-invariant unobservables (z _i and z _t), but latent attributes could change over time and their significance could also vary (Noel and Nyhan, 2011); that is, unobserved time-varying confounding, $z_{it} = a_t\ast z_i$, is still left in the error term.

Network scholars have considered the factor approach as an identification strategy. For instance, when Noel and Nyhan (2011) discuss unobserved homophily with time-varying impact, they suggest such homophily is ”likely to be a factor” (p. 212); Van der Weele (2011) explicitly specifies unobserved homophily as a latent factor U in sensitivity analysis of network effects; and He and Hoff (2019) develop “simple and scalable linear regression and latent factor models” (p. 2) as an alternative to SAOMs and model unobserved nodal attributes as latent factors. But the behavior equation of their models estimates exogenous rather than endogenous network influence. In this research, we advance the idea of using multiple factors as an identification strategy and extends the existing spatial-temporal regressions to be a multifactor model.

We develop a Markov Chain Monte Carlo (MCMC) algorithm to estimate the proposed MLST-MF model and adopt Bayesian shrinkage to determine the number of factors. Monte Carlo studies find (1) Bayesian shrinkage is able to correctly detect the number of latent factors even with uninformative networks and strong confounding; (2) compared to two-way fixed effects, the multifactor term significantly improves identification of network interdependence; (3) the varying-interdependence specification not only recovers temporal heterogeneity of network influence but also detects invariant or null interdependence; and (4) falsely assuming interdependence to be time-invariant results in highly misleading inferences. We then apply the proposed method to two examples in international relations (IR) to illustrate how it can be used in empirical studies. An R package bpNet is developed to assist researchers to implement this method.

The proposed method is squarely situated in the literature of network effect analysis. We take the network selection process into consideration but treat it as a nuisance, though studies on network dynamics have formed a large and growing literature. Our approach can be viewed as a simple alternative to a more comprehensive solution that coevolution modeling may offer. Our method models endogenous network effects and handles latent confounders, two important tasks that are largely neglected in the coevolution literature. Also, unlike coevolution models, our model does not rely on the strong modeling assumptions about the selection process or impose restrictions on network ties (Windzio, 2021). Therefore, the MLTS-MF model is robust to mis-specification and flexible to accommodate various types of networks.

Our method also has several limitations. First, it requires a large number of units or time periods to precisely estimate the latent factors or factor loadings. Second, the estimates of latent factors are not substantively interpretable and no inference can be drawn about the important missing attributes. Third, we only consider continuous and monodic outcomes. Future research may extend the model to a generalized linear model or further deal with methodological issues in modeling dyadic outcomes. Finally, this research conducts identification in the regression framework. The simultaneous nature of network interdependence does not fit the conventional causal inference approaches, and the framework of chain graphs may be a future direction to pursue identification in a strict sense of causality (Tchetgen et al., 2021; Ogburn et al., 2020).

1. The model

Suppose there are N units and each unit i is repeatedly measured at time t = 1,  2,  …,  T. The MLST-MF model is specified as below:

(1)$$y_{it} = \rho_t{\bf w}_{it}{\bf y}_t + {\bf w}_{it}{\bf X}_t {\boldsymbol \beta}_1 + \gamma y_{i, t-1} + {\bf x}_{it}{\boldsymbol \beta}_2 + \nu_i + \psi_t + {\boldsymbol \zeta}_{i} {\bf f}_t + \epsilon_{it},\; $$

(2)$$\rho_{t} = \alpha_0 + \kappa\rho_{t-1} + {\bf Z}_t {\boldsymbol \alpha} + \eta_{t},\; $$

where y _it is a monodic outcome of ego i at time t. In Equation (1), the spatial autoregressive coefficient ρ_t is the quantity of interest and measures the strength of outcome interdependence in network W_t. The vector w_it is the ith row of W_t, and W_t is a row-standardized (N × N) adjacency matrix with diagonal elements equal to zero. Here we impose no restrictions on the network: It can consist of well-separated groups or fuzzy clusters, and its ties can be binary or continuous, directed or undirected, changing or constant. The vector y_t = {y _1t,  …,  y _jt,  …,  y _Nt} includes outcomes of all units at t. Because the outcome variable appears on both sides of Equation (1), ρ_t is determined by simultaneous moves made by all agents in the network. We allow ρ_t to vary over time and assume an autoregressive process with order 1 as in Equation (2). The process also includes group-level covariates Z_t to further model ρ_t's variation. Other options of the ρ_t process are also possible, such as a random walk process for local smoothing (ρ_t = ρ_t−1 + Z_tα + η_t) or a varying-slope specification (ρ_t = α₀ + Z_tα + η_t). Based on our tests, these options give very similar estimates of ρ_t. Technically, Z_t can be any unit-invariant variables. When Z_t includes structural features of W_t (e.g., network density, connectivity, reciprocity, etc.), network interdependence varies with nodal behavior at the micro level and structural changes at the macro level. The second network term w_itX_tβ₁ is to estimate exogenous network effects, where X_t = (x_1t,  …,  x_Nt) may include alters’ attributes or/and previous outcomes.²

The model has three components that help identify ρ_t. First, we follow Christakisa and Fowler and include ego's lagged outcome y _i,t−1. The term γy _i,t−1 can capture confounders whose impact is mediated by the previous outcome, blocking the backdoor path, Path in Figure 1, y _jt ← y _j,t−1 → w _ijt ← y _i,t−1 → y _it. The term also helps to correct for serial correlation. Second, x_it β₂ includes observed nodal attributes and contextual factors to block y _jt ← (x _j,   x _j,t−p,  S _1t) → w _ijt ← (x _i,   x _i,t−p,  S _2t) → y _it, Path and Path in Figure 1. Third and most importantly, the multifactor term ζ_i f_t is to approximate latent confounders, where f_t is a (r × 1) vector of factors, and ζ_i is a (1 × r) vector of factor loadings. The number of factors, r, will be determined by Bayesian shrinkage, which we will discuss later. The multifactor term is to block $y_{jt}\leftarrow ( \alpha _j,\; \, \ \alpha _{j, t-p},\; \, f_{1t}) \rightarrow w_{ijt} \leftarrow ( \alpha _i,\; \, \ \alpha _{i, t-p},\; \, f_{2t}) \rightarrow y_{it}$, Path and Path in Figure 1. We restrict the first factor and the second factor loading to be 1, i.e., f_t = (1,  ψ_t,  …,  f _rt) ′ and ζ_i = (ν_i,  1,  …,  ζ_ri), and the multifactor term includes two-way fixed effects. In Equation (2), we write ν_i and ψ_t as separate terms from the multifactor term to emphasize that the model also considers time- and unit-invariant confounders.

Fig. 1. Identifying network interdependence: backdoor paths and sources of confounding.

We use factor analysis to reduce bias based on the following assumption: for each node i = 1,  2,  …,  N in a longitudinal network, there exists a set of unobserved time-varying confounders U_i = (u _i1,  u _i2, · · · ,  u _iT), such that the stacked (n × T) matrix U = (U₁,  …,  U_n) can be approximated by two lower-rank matrices (r ≪ min{N,  T}), i.e., U = L′F in which F = (f₁,  …,  f_T) is a (r × T) matrix of factors and L = (ζ′₁,  …,  ζ′_n) is a (r × N) matrix of factor loadings (Pang et al., 2021).

To determine the appropriate number of factors is crucial to ensure sufficient bias-reduction but no overfitting. There are several methods for factor selection, including model comparison (Pang, 2014; Bai and Li, 2021), cross-validation (Xu, 2017), Bayesian shrinkage techniques (Kyung et al., 2010; Pang et al., 2021), and so on. We adopt the shrinkage approach for its computational efficiency. Because the standard Bayesian LASSO prior is not directly applicable to the varying parameters in the model, we use the hierarchical shrinkage method proposed by Bitto and Frühwirth-Schnatter (2019). The LASSO-like method works as the following. We re-parameterize factor loadings ζ_i as ${\boldsymbol \zeta }_i = {\boldsymbol \omega }_{{\boldsymbol \zeta }} \cdot \tilde {{\boldsymbol \zeta }}_i$ for all i, where ${\boldsymbol \omega }_{{\boldsymbol \zeta }} = ( \omega _{\zeta _1},\; \, \ldots ,\; \, \omega _{\zeta _{r}})$. Assuming the variance-covariance matrix of the factor loadings ζ as ${\bf \Omega }_{\zeta } = \rm {Diag}( \omega _{\zeta _1}^{2},\; \, \ldots ,\; \, \omega _{\zeta _{r}}^{2})$, the variance-covariance matrix of $\tilde {{\boldsymbol \zeta }}_i$ is an identity matrix. Then we apply the LASSO prior on each element in ω_ζ. When $\omega _{{\boldsymbol \zeta }_{m}}$ is shrunk to zero, the mth factor loading vector ${\boldsymbol \zeta }_{m} = \omega _{{\boldsymbol \zeta }_m} \cdot \tilde {{\boldsymbol \zeta }}_{m}$ is also shrunk to zero, and its associated mth factor f_m = {f _m1,  …,  f _mT} is excluded from the model.

Bayesian shrinkage lacks sparseness and will never generate estimates exactly equal to zero. Redundant factors are only virtually excluded with negligible impact on the outcome. Which factors are included or virtually excluded can be easily observed through the estimated posteriors of ω's. Because ω² = ( ± ω) ², the posterior distribution of ω should be clearly bimodal when its associated factor successfully escapes shrinkage. But the posterior distribution of ω will be unimodal and centered at zero when its associated factor is virtually excluded from the model. When the data are undecisive about whether to include a factor, the posterior of the associated ω will be a mix of a spike around zero and a bimodal distribution.

Because factor analysis is a dimension-reduction technique but not an estimation of latent nodal attributes or common exposures, Bayesian shrinkage simply chooses the number of factors to best summarize the residual covariance. Neither the shrunk factors nor the factors that successfully escape shrinkage can be interpreted as approximations of specific attributes or common shocks, unless we have strong theories about both the selection and influence processes. Therefore, the factor and factor loadings should be treated as nuisance parameters that are used to correct for bias but not as estimates for making inferences about the network evolution.

Finally, we assume that the errors in the individual- and group-level regressions are independently, identically, and normally distributed: $\epsilon _{it} \overset {iid}\sim {\cal N}( 0,\; \, \sigma _e^{2})$ and $\eta _t \overset {iid}\sim {\cal N}( 0,\; \, \sigma _{\eta }^{2})$. To complete the Bayesian model specification, we assign priors to parameters in the model. The functional forms of priors are reported in the online Appendix A.1.1. We develop a MCMC algorithm to estimate parameters and select factors. Time-varying network interdependence, ρ_t, is updated with the Metropolis–Hastings algorithm with tailored proposal densities, and the rest of the parameters are updated with the Gibbs sampler. We report the formal expression of the posterior distribution and the MCMC sampling procedure in Appendix A.1, where we also discuss the stationarity issue.

2. Monte Carlo studies

We investigate the performance of the MLST-MF model on simulated data and test how the two major extensions—the varying-ρ_t specification and the multifactor error structure—improve identification of network interdependence. This section reports three Monte Carlo studies. The first exercise is designed to be a tough test considering a network with little variation. In the second study, we set network interdependence to be exactly zero to see whether the proposed varying-ρ_t model is able to detect invariant and null interdependence. The third study is intended to further reveal the consequences when the constant interdependence assumption is violated, which demonstrates the importance of allowing the spatial autoregressive coefficient to vary.

In each study, we compare four models: M1 is MLST-MF including the true number of factors, M2 is MLST-MF including ten factors and using Bayesian shrinkage, M3 is a fixed-effect MLST with no factors, and M4 has the same specification as M1 except that network interdependence is set to be constant. Table 1 summarizes features of data generating processes (DGPs) and model components. In Study I and Study II, the true network interdependence is set to be ρ_t = 0.3 + 0.3ρ_t−1 − 0.2 Z _t + η_t, $\eta _{t} \overset {iid}{\sim } {\cal N}( 0,\; \, 0.36)$. In these studies, we also test how model performance varies when the tie-homophily correlation gets stronger. We discuss estimates of ρ_t (or ρ) and results of factor selection in the main text. Posteriors of other parameters are reported in the online Appendix.

Table 1. DGPs and model specifications in monte carlo studies

2.1 Study I

We generate a time-invariant network W that consists of 20 separate groups each with ten members. When two units i and j belong to the same group, their ties are w _ij = w _ji = 1; otherwise, w _ij = w _ji = 0. This network lacks variation and makes it particularly difficult to separate network effects from other sources of correlated outcomes (Lin, 2010). To make identification even more challenging, we generate confounding by sampling the factor loadings ζ from a multivariate normal distribution ${\cal N}( {\bf 0},\; \, \kappa _w {\bf \Sigma })$. The off-diagonal elements of Σ are the same as in W and the diagonal elements of κ_wΣ are set to be 1. The parameter κ_w controls the strength of confounding: there is a tie between unit i and unit j (w _ijt = 1) if cor(ζ_i,  ζ_j) ≥ κ_w; otherwise, w _ijt = 0. We set κ_w = 0.3,   0.6,   0.9, respectively, to generate three datasets. The DGP of the outcome y _it is reported as Equation (A24) in the online Appendix.

Figure 2 depicts the estimated trajectory of ρ_t (the dashed line) with the 95% credibility interval (CI, the gray band) based on the four models, plotted against the true time series of ρ_t (the black solid line). Figure 2 (a) shows M1 performs very well in recovering the true ρ_t's at all levels of confounding. Only when the tie-homophily correlation is as strong as κ_w = 0.9 is ρ_t sometimes slightly over-estimated. M1 includes two factors a priori, but in reality we do not know the true number of factors. M2 relies on data and algorithm to determine the number of factors. Figure 2 (b) shows that M2 performs almost as well as M1 when confounding is weak or modest, and overestimates ρ_t more often than M1 when κ_w = 0.9. Figure 3 (a) reports the posterior distributions of ω's when κ_w = 0.9. There are eight spike posteriors centered at zero and two bimodal ones, clearly suggesting two factors out of the ten initial factors. Bayesian shrinkage also correctly detects the number of factors when κ_w = 0.3,  0.6. Due to space limitations, we do not report the posteriors of ω's in these simulations.

Fig. 2. Posteriors of network Interdependence (Study I). Note: In each graph, the solid line indicates the true time series of ρ_t, the dashed line is the posterior mean over time, and the gray area is the 95% Bayesian CI.

Fig. 3. Factor selection. Note: The bimodal posteriors indicate that their associated factors should be included, while the unimodal ones indicate that their associated factors are virtually excluded from the model.

Figures 2(c) shows that M3 performs much worse than M1 and M2. The model simply uses fixed effects to correct for bias and overestimates ρ_t even when confounding is weak (κ_w = 0.3). Estimates move farther and farther away from the true values when the tie-homophily correlation gets stronger. Finally, M4 is the model which mis-specifies the varying network interdependence as a constant ρ. It neither estimates the time heterogeneity of interdependence nor produces a sensible estimate of the time-average of ρ_t, as shown in Figure 2(d). In fact, the posterior of ρ could not converge, indicating mis-specification.

This exercise demonstrates that the proposed model performs well even when identification is very challenging with an uninformative network and strong confounding. The Bayesian shrinkage method is able to correctly detect the number of latent factors. The multifactor approach significantly improves identification compared to the two-way fixed effects and effectively reduces bias caused by time-varying confounders.

2.2 Study II

Study I shows that the ρ_t-MLST-MF model performs well when it is correctly specified. But how does the model perform when the true interdependence is invariant and even exactly zero? In Study II, we answer the question by setting $\rho _t = 0,\; \, \ \forall t$ and re-running the simulations in Study I.

Figure 4 (d) displays that the correctly specified model, M4, performs well—the estimated posterior of ρ contains zero in all the three simulations. The other three varying-ρ_t models are all misspecified, but they perform differently. The posteriors based on M1 and M2 correctly suggest that network interdependence barely varies and is likely to be zero. The estimated trajectories in Figure 4 (a) and (b) vary within a very narrow band centered at the level close to zero with the 95% CI containing zero most of the time. The only exception is M2 when κ_w = 0.9, but the bias is rather small. Figure 3 (b) reports the results of factor selection when κ_w = 0.9, and Bayesian shrinkage again correctly detects the number of latent factors to be 2. Figure 4 (c) shows that M3 mistakenly suggests positive and significant network interdependence, and bias gets more serious when the homophily-tie correlation gets stronger. Furthermore, the posterior mean of ρ_t in M3 varies much more widely than in M1 and M2, falsely suggesting salient temporal heterogeneity of interdependence.

Fig. 4. Posteriors of network interdependence (Study II).

These results show that MLST-MF with varying ρ_t is a robust specification. The robustness comes from the multifactor specification because the fixed-effect version of the model performs poorly when mis-specified.

2.3 Study III

In Study III, we further investigate the consequences when the time-invariant interdependence assumption is violated. We make identification relatively easy in order to focus on the impact of the specification of interdependence on model performance. We generate data following the exact DGP as in Study I but using an empirical network. The network is constructed with ICEWS dyadic event data (He and Hoff, 2019) in 50 time periods, consisting of the 50 most active countries. It varies frequently over time and has continuous and asymmetric ties. These features help identify and estimate network interdependence. Also, there are no latent confounders in this DGP. For the purpose of this study, we focus on comparing M2 and M4.

As shown in Figure 5(a), the ρ_t-MLST-MF model with Bayesian shrinkage on teninitial factors recovers the true time series of ρ_t. The CIs are wider than those in the previous studies only because the sample size is smaller. The ρ-MLST-MF model includes the true number of factors a priori, but generates a point estimate of ρ that is much larger than the average level of the true ρ_t's. As shown in Figure 5(b)-(i), the posterior lies far above 0.3. Figure 5(b)-(ii) depicts the posterior kernel of ρ. Unlike in Study I, the Markov Chain does not clearly display non-convergence or poor mixing. No sign of numeric problems makes it difficult to detect mis-specification, and we are likely to make highly biased inference about network interdependence.

Fig. 5. Posterior of network interdependence.

This exercise highlights the necessity and importance of relaxing the assumption of invariant interdependence because a violation of this assumption incurs serious bias even when there are no latent confounders.

The three simulation studies together reveal several important findings. First, the multifactor analysis is an effective bias-reduction method and performs well in hard cases when the network is not informative, confounding is strong, and the true network interdependence is constant or even zero. Second, falsely assuming constant interdependence has serious consequences. It not only neglects temporal heterogeneity of ρ_t but also mis-estimates the central level of ρ_t. Third, Bayesian shrinkage is able to detect the true number of factors. Leaving the redundant factors in the model slightly affects identification when the network is not informative. Fourth, network types and data structure notably affect identification. Identification is particularly difficult with time-invariant, binary, and symmetric networks (Shalizi and Thomas, 2011). The number of units or time points also matter, because latent factors or factor loadings are more precisely estimated with a larger spatial or time dimension.

3. Empirical applications

We apply MLST-MF to two empirical IR studies to demonstrate the implementation of the proposed method.

3.1 Migration and terrorism

Bove and Bohmelt (2016) (henceforth, B&B) investigate whether the global immigration network leads to clustered terrorist attacks, drawing on data covering 145 countries between 1970 and 2000. The outcome variable is the number of terrorist attacks (on a logarithm scale) in the territory of a country in a year. The network tie is defined as the number of immigrants (refugees excluded) to country i from country j in year t. The network data vary considerably at the turn of each decade but very little in the rest of the sample years. To study the effect of a longitudinal network, B&B still assume invariant network influence. The authors discuss the possibility that correlated terrorist attacks could be “a mere product of a clustering in similar [state] characteristics (p. 578).” They adopt two strategies to avoid spurious network interdependence. The first strategy is to include country- and year-dummies in SAR, which finds the exogenous network effect is 0.08 (Model 3 in the original paper). The second strategy relies on the m-STAR model (Hays et al., 2010) and uses geographical adjacency to approximate omitted homophily. The endogenous network effect is estimated as 0.07 (Model 5 in the original paper). Both estimates are statistically significant and support their theory.

We re-analyze the data with three different specifications, including the full model (ρ_t-MLST-MF), the fixed-effect model (ρ_t-MLST-FE), and the constant-interdependence model (ρ-MLST-MF). These model specifications include the same set of covariates as in the original study, but they also add several components that are not in B&B, including a lagged outcome, multiple latent factors, an exogenous network-effect term of $\beta ( {\bf W}_t\ast Domestic \ Conflict)$, and the network density as a group-level predictor. The definitions and descriptive statistics of the variables are reported in Table A1 in the online Appendix.³

When estimating ρ_t-MLST-MF and ρ-MLST-MF, Bayesian shrinkage suggests 1–4 latent factors out of the ten initial factors (Figure A8). Compared to the simulation examples, empirical data are more ambiguous about how many factors to be included mainly because factors could be highly correlated in real data analyses. Figure 6 reports the estimated network interdependence based on the three models. All models suggest positive network interdependence, and these estimates are statistically reliable most of the time, except for results based on ρ_t-MLST-MF in the years after 1992. As for the general strength of interdependence, ρ_t-MLST-FE estimates the time-average of ρ_t as 0.079, almost the same as the estimate $\hat {\rho } = 0.08$ based on SAR with country and year dummies in the original study. But the estimate is much smaller according to the constant-ρ model (ρ-MLST-MF), and the posterior mean of ρ_t is 0.066, close to what the m-STAR model finds ($\hat {\rho } = 0.07$) in B&B. The full model, ρ_t-MLST-MF, estimates the time series of ρ_t that varies around the level of 0.060. Because omitted homophily and common shocks tend to generate upward biases (Van der Weele, 2011), the smallest estimate based on ρ_t-MLST-MF suggests that the full model corrects for bias most sufficiently among all the models considered.

Fig. 6. Estimated network interdependence (migration and terrorism).

The finding about the temporal variation of interdependence is interesting: ρ_t changes non-monotonically. Figures 6 and 6(b) show that the two varying-ρ_t models agree with each other and suggest that the trajectory of ρ_t is clearly inverse V-shaped—the effect increases in the first decade and then decreases over time. The two models disagree on the level of the ρ_t trajectory: the curve of posterior mean in Figure 6(b) universally dominates the curve in Figure 6 (a).

Our analysis supports the key hypothesis in B&B: the global immigration network transmits the risk of terrorist attacks. But the proposed model finds the effect is much smaller and changes nonlinearly over time. The posteriors of the coefficients in the individual- and group-level regressions are reported in Figure A9 in Appendix A.3.1. They mostly confirm the estimates reported in B&B.

3.2 GATT/WTO and free trade

The second example is Chaudoin et al. (2015) (henceforth, CMP) who propose a varying-ρ_t SAR model and apply it to estimate interdependence of trade policies of developing countries in the GATT/WTO network. Scholars have been debating about whether the wave of trade liberalization in the developing world since the 1970s is primarily attributed to changes in domestic politics (e.g., democracization) or mainly explained by the international system (e.g., the global trade institution of GATT/WTO). From the methodological perspective, this debate urges a separation of the effect of homohily (i.e., joining GATT/WTO and lowering trade barriers are both consequences of a country's democratization) from network interdependence (i.e., the mutual effect of GATT/WTO members’ trade policy choices).

CMP focus primarily on the effect of democratization. Although they are not substantively interested in network interdependence, CMP relax the constant interdependence assumption and include two-way fixed effects in their extended SAR model. Based on a sample of 85 developing countries between 1970 and 2008, CMP find democratization significantly reduces tariff rates and drives developing countries to move together to free trade. In contrast, tariff rates have been negatively interdependent among GATT/WTO member states except in a few years in the 1980s, and the strength of policy interdependence varies widely during the sample years (Table 4 in the original study, and Figure A10(a) in our online Appendix).

For the purpose of this study, our reanalysis focuses on identifying and estimating the GATT/WTO network influence. The network is binary and symmetric: w _ijt = w _jit = 1 if countries i and j are both GATT/WTO member states in year t; w _ijt = w _jit = 0, otherwise. The network has several special features, including (1) it consists of an expanding clique of GATT/WTO member states and a decreasing number of non-member states as isolated nodes; (2) network links never dissolve; and (3) when a non-member state forms network links, it does not choose individual “neighbors” but joins the clique, forming links with all existing member states. These tie-formation features make SAOMs inapplicable because SAOMs require a delicate balance between the stability and change of the network. However, this type of network is common in IR studies because IR scholars are often interested in the influence of multilateral insinuations on the diffusion of national policies.

In the original study the outcome variable is the tariff rate of a country in a year. The variable is bounded to be non-negative and the data are highly skewed, which violates the normal error assumption. Our initial analysis suggests 30 factors, which is a sign of mis-specification.⁴ Therefore, we take a logarithm of tariff rates. Besides transforming the dependent variable and using multiple factors, our specification is different from the CMP model in several other respects: we add a lagged dependent variable, an exogenous network effect term,⁵ and a state-equation of ρ_t = ρ₀ + αZ _t + η_t, where Z _t is the yearly density of the GATT/WTO network. We include the same set of covariates as in CMP, and their definitions and descriptive statistics are reported in Table A2 and Table A3 in the online Appendix.⁶ Bayesian shrinkage suggests 5–7 latent factors, among which 5 should be certainly included, as shown in Figures A12 and A13.

Figure 7 depicts estimated posteriors of ρ_t's based on three models, ρ_t-MLST-MF, ρ-MLST-FE, and ρ-MLST-MF. They all suggest that GATT/WTO drives national trade policies of developing countries to diverge rather than converge. The negative network interdependence is puzzling, but confirms the finding in the original study. According to ρ-MLST-MF, the posterior mean of ρ is − 0.196 and the 95% CI do not contain zero, as shown in Figure 7(c). The two varying-ρ_t models also find negative interdependence, but the estimated trajectories are notably different from each other. The multifactor model finds ρ_t slightly decreases (i.e., the negative interdependence gets stronger) over time, whereas the fixed-effect model suggests ρ_t is in a steady increase (i.e., the negative interdependence gets weaker over time). Also, ρ_t changes much slower according to ρ_t-MLST-MF than ρ_t-MLST-FE. Figure 7(a) shows that posterior mean of ρ_t varies from − 0.146 in 1972 to − 0.195 in 2008, and Figure 7(b) indicates the correspondent values are − 0.273 and − 0.085. These findings suggest that the GATT/WTO network effect on tariff rates is negative and relatively stable over time.

Fig. 7. Estimated network interdependence (GATT/WTO and free trade).

The ρ_t-MLST-FE model generates very different estimates of the other coefficients from those based on the multifactor models, as shown in Figure A14 in the online Appendix. All three models estimate the effect of democratization (Regime) to be negative and significant. This echos the finding in CMP, showing that democratization drives developing countries to move together to free trade despite the dispersing effect of the GATT/WTO institutional network. This empirical finding highlights “the need for improvements in our theoretical arguments for how WTO membership affects members’ tariffs directly and through contagion across members (CMP, p. 299).”

4. Conclusion

This paper proposes a multilevel spatio-temporal model with a multifactor error structure to estimate temporally heterogeneous network interdependence. Studies on simulated and empirical data demonstrate that the multifactor method with Bayesian shrinkage is an effective bias-reduction strategy. Also, a violation of the constant interdependence assumption leads to biased and misleading inferences about network interdependence. The Bayesian model is highly flexible and can accommodate various types of networks. Besides endogenous network effects, MLST-MF can also be applied to identify and estimate exogenous network effects since exogenous and endogenous network effects are subject to the same identification problem. When the researcher is only interested in exogenous network effects, she can apply a simplified version of the proposed model without the spatial lag term.

The MLST-MF model provides a simple alternative to a more comprehensive solution that may directly model the interaction of the selection and influence processes. The limitations and restrictions of the existing co-evolution models make alternatives such as MLST-MF needed and useful. Table 2 summarizes the features of MLST-MF compared to coevolution SAOMs. MLST-MF has its advantages and disadvantages compared to SAOMs, and is more appropriate than the SAOM model in situations where (1) network effects are the primary quantities of interest; (2) selection is not substantively interesting but mainly concerned as a source of bias; (3) it cannot be well justified that there are no important confounders left in the error term. But the researcher should choose coevolution models over our model when the mutual effects of selection and influence are substantively interesting and when there is no need to worry about latent homophily and common exposures. In the online Appendix, we report simulated and empirical examples to compare MLST-MF and SAOM and illustrate how MLST-MF can be applied as an alternative to SAOM. Future research may explore the idea of selection-correction co-evolution models (e.g., Han et al., 2021), probably by incorporating multifactor terms in selection and behavior equations to reduce bias from mis-specification of the selection process and latent confounders in the influence analysis.

Table 2. Comparison of the MLST-MF model and SAOM-based co-evolution models.