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.Footnote 1 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, Reference Elkins and Simmons2005), and strategic interaction and optimizing behavior of networked actors are deemed as the sources of outcome interdependence (Manski, Reference Manski1993; Franzese and Hays, Reference Franzese and Hays2008b). Network interdependence is substantively interesting because it generates social multipliers and has unique implications for policy interventions (Burke, Reference Burke2016).
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., Reference Bramoulle, Kranton and D'Amours2014). 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, Reference Xu2017; Samartsidis et al., Reference Samartsidis, Seaman, Montagna, Charlett, Hickman and Angelis2020; Pang et al., Reference Pang, Liu and Xu2021). 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, Reference Manski1993) and the inverse Galton's problem (Gleditsch and Ward, Reference Gleditsch and Ward2006). The challenge can be illustrated with the directed acyclic graph approach (Pearl, Reference Pearl2009): 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, Reference Shalizi and Thomas2011; Elwert and Winship, Reference Elwert and Winship2014). 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., Reference Han, Hsieh and Kox2021), and endogenous network influence is seldom the quantity of interest in such models. More importantly, coevolution models are still subject to latent confounders (Robins, Reference Robins2015; Ragan et al., Reference Ragan, Osgood, Ramirez, Moody and Gest2019; He and Hoff, Reference He and Hoff2019). Furthermore, they impose restrictions on network ties and are applicable primarily to networks built on binary ties (Snijders, Reference Snijders2001, Reference Snijders2005; Steglich et al., Reference Steglich, Snijders and Pearson2010). 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, Reference Christakisa and Fowler2008, Reference Christakisa and Fowler2013). This argument has been widely challenged and their model also suffers from other estimation and inference problems (Robins, Reference Robins2015). 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, Reference Noel and Nyhan2011); 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 (Reference Noel and Nyhan2011) discuss unobserved homophily with time-varying impact, they suggest such homophily is ”likely to be a factor” (p. 212); Van der Weele (Reference Van der Weele2011) explicitly specifies unobserved homophily as a latent factor U in sensitivity analysis of network effects; and He and Hoff (Reference He and Hoff2019) 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, Reference Windzio2021). 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., Reference Tchetgen, Fulcher and Shpitser2021; Ogburn et al., Reference Ogburn, Shpitser and Lee2020).
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:
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 Wt. The vector wit is the ith row of Wt, and Wt 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 yt = {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 Zt 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 + Ztα + ηt) or a varying-slope specification (ρt = α0 + Ztα + ηt). Based on our tests, these options give very similar estimates of ρt. Technically, Zt can be any unit-invariant variables. When Zt includes structural features of Wt (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 witXtβ1 is to estimate exogenous network effects, where Xt = (x1t, …, xNt) may include alters’ attributes or/and previous outcomes.Footnote 2
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, xit β2 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 ft is to approximate latent confounders, where ft 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., ft = (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.
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 Ui = (u i1, u i2, · · · , u iT), such that the stacked (n × T) matrix U = (U1, …, Un) can be approximated by two lower-rank matrices (r ≪ min{N, T}), i.e., U = L′F in which F = (f1, …, fT) is a (r × T) matrix of factors and L = (ζ′1, …, ζ′n) is a (r × N) matrix of factor loadings (Pang et al., Reference Pang, Liu and Xu2021).
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, Reference Pang2014; Bai and Li, Reference Bai and Li2021), cross-validation (Xu, Reference Xu2017), Bayesian shrinkage techniques (Kyung et al., Reference Kyung, Gill, Ghosh and Casella2010; Pang et al., Reference Pang, Liu and Xu2021), 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 (Reference Bitto and Frühwirth-Schnatter2019). 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 fm = {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 ω2 = ( ± ω) 2, 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.
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, Reference Lin2010). 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.
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.
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.
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, Reference He and Hoff2019) 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.
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, Reference Shalizi and Thomas2011). 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 (Reference Bove and Bohmelt2016) (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., Reference Hays, Kachi and Franzese2010) 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.Footnote 3
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, Reference Van der Weele2011), the smallest estimate based on ρt-MLST-MF suggests that the full model corrects for bias most sufficiently among all the models considered.
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. (Reference Chaudoin, Milner and Pang2015) (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.Footnote 4 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,Footnote 5 and a state-equation of ρt = ρ0 + α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.Footnote 6 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.
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., Reference Han, Hsieh and Kox2021), 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.
Comparison of the MLST-MF model and SAOM-based co-evolution models.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/psrm.2022.47. and Replication https://doi.org/10.7910/DVN/B5RVWB
Acknowledgments
These authors contributed equally to this work. We thank Nick Beauchamp, Mathew Blackwell, Chong Chen, Bruce A. Desmarais, Naoki Egami, Jeff Gill, In Song Kim, Teppei Yamamoto, Luwei Ying, and Danlin Yu, and all participants at the Panel on Spatial Analysis at the virtual annual conference of Society for Political Methodology in Toronto on 15 July 2020, the Networks Panel at the American Political Science Association annual meeting on 31 August 2019, and the MIT Political Methodology Speaker Series on 29 October 2019. We are especially grateful to Associate Editor Daniel Stegmueller and two anonymous reviewers for their incredibly helpful comments. Xun Pang acknowledges support from the National Social Science Foundation of China (Grant No.17ZDA110).
