5.1. GNN Models

GNNs are a generalization of neural networks to graph data (Labonne Reference Labonne2023). GNNs encode representations of nodes and edges that depend on the graph structural properties, as well as node and edge features. A GNN functions through neural message passing which passes information between nodes and is updated through the neural networks. One of the most versatile forms of GNNs is the Graph Convolutional Network (GCN) architecture (Kipf and Welling Reference Kipf and Welling2017). This is an adaptation of Convolutional Neural Networks (CNNs) usually used for deep learning about images and texts to graphs (Zhou et al. Reference Zhou2020). A GCN is specifically defined with a layer-wise propagation rule, which usually is a non-linear function that updates node and edge features by aggregating the features of their neighbors. GCNs efficiently leverage the graph structure by performing convolutions, thereby enabling effective feature extraction for tasks like node classification and embedding, edge embedding and prediction, and graph embedding and classification. In this article, I first use a deep GNN model which includes multiple layers in the model propagation as shown in the generic example of Figure 4.

Figure 4

GCN model (Kipf and Welling Reference Kipf and Welling2017).

Generally, there are three types of layers: input, hidden, and output. These layers are chosen depending on the task at hand. In this article, I choose as the basic GNN layer for purpose of graph classification what is known as GAT to incorporate not only node features but also (as explained before) various edge features. GATs make use of attention processes to give neighboring nodes varying weights depending on their relative relevance in the graph (Zhou et al. Reference Zhou2020), with the weights learned during the training phase. This enables the model to zero in on the most relevant information in the graph. GAT thus dynamically enhances the ability of the model to capture highly complex interactions in graph-structured data. I specifically use the EdgeGATConv layer (see Supplementary Material) as the base for constructing three different architectures geared, respectively, toward (a) graph classification using topological labels, (b) conditioned extraction of latent features, and (c) graph generation.

5.2. Graph Classification

The topological information stored in the FPCA is discretized into four categories which are then used to classify the graphs into four topological classes through supervised deep machine learning.Footnote ⁷ The list of graphs is split into training and testing sets of graphs.Footnote ⁸ The goal is to calculate the accuracy of predicting graph distribution after training the GNN classification model.

Having trained and saved the model, we can use the saved model to examine, for example, how the membership probabilities in topological classes change from graph to graph, that is, how the topological structure of the international system of states evolves from year to year.

Figure 5 displays the yearly proportions of topological classes. The colors indicate the different topological classes from the data. The boundaries between the class probabilities change over the years, which means that the topological structure of the international system varies over the years, with most variation occurring between the mid-1960s and the late 1980s, the latter corresponding to the end of the cold war, with Class 0 being predominant by far. During these twenty years or so the international system saw many changes in addition to the re-ignition of the cold war. We then see important shifts in the boundaries, with Class 3 gaining in membership whereas the others shrinking, if moving toward stabilization, during the 1995–2010 period, with Class 2 being the dominant. The few years in the aftermath of the cold war, which witnessed much instability in world order, are manifested in the sharp variation of the probabilities for the four classes, with Class 3 being dominant during that period.

Figure 5

Proportions of clusters over the years.

So far, these plots have been done using the predictive power of a graph classification model. However, the basic layers of the model utilize, as explained before, the attention mechanisms, the weights of which can be learned. These weights are edges weights (not nodes) and hence we can use them to extract some useful information concerning how edge importance is propagated in the neighborhoods (which is what the attention mechanism focuses on). I first look at the global picture by focusing on the whole set of graphs (that is, from 1946 to 2010). Figure 6 shows how the maximum attention weights vary in conjunction with persistence entropy and total MIDs in the system (all three have been spline-smoothed for comparison clarity purposes).

Figure 6

Max attention weights, persistence entropy, and MIDs.

The time trends of the three quantities are quite similar in general. This signifies that the attention mechanisms are also able to capture the dynamics in terms of MIDs, much like persistence entropy does. Persistence entropy is a measure of topological disorder, total MIDs per year can be understood as a measure of political disorder in world politics, the evolution of both is roughly speaking picked up by the maximum attention weights of the edges (relations between states). The upward behavior of the total MIDs curve is not captured by either the persistence entropy or the attention weights.Footnote ⁹

Second, we can also zoom in by finding which edges have the top k attention weights in any single year (graph). Plotting these, we obtain Figure 7 which displays, at the top, the subgraphs with the eight highest attention weights for 1985 and 2010. The number displayed in the edges are the normalized attention weights between the respective nodes (states). Not only does the structure of the subgraph change but also the nodes (states) composing the subgraphs. The lack of an edge between any two nodes means that there isn’t any meaningful attention propagation between them. Remember that these attention mechanisms provide information on which edges are important in their neighborhoods, and how information is passed through the neighborhood. Edge-level attention focuses on assigning weights to edges rather than nodes. This is particularly, useful in scenarios where the relationships between nodes (represented by edges) carry different levels of importance. In the empirical illustration, much information consists of features of the edges in a graph, and hence attention weights allow us to learn the importance of these features in the model, both locally and globally. The bottom two plots represent the subgraphs for the same eight states with the normalized number of shared IOs as a weight (red) in 1985 and 2010. These IO-sharing membership weights are not learned; they are computed from the sociomatrix. The attention weights are optimally learned, which means that they capture the dynamic relations as the neural networks are fired throughout training, that is, they reflect the dynamic and highly complex short and long “range” correlations in the data across the neural networks. While not done in here, attention weights can be used for classification as well as prediction purposes. They can also be used to draw heatmaps of the various topological classes as a function of time, thereby providing useful information on the evolution of the system of graphs.

Figure 7

Subgraphs with top eight attention and IOs-sharing weights.

All results obtained in this section, used the TDA-generated topological measures as “classification labels”. In the next section, I construct a model which uses the topological effects as a constraint that conditions the space of latent factors. This is an extension of the idea of using topology effects as moderators in the Bayesian model as explained earlier.

5.3. GNN-Extracted Latent Features

Analysis of latent factors is used in statistical analysis to reduce dimensionality, improve model accuracy, and uncover relationships not immediately observed, using tools such as Confirmatory Factor Analysis (CFA) and SEM. Latent factors are important for understanding causality by controlling for unobserved con-founders, thereby leading to more precise conclusions. Because GNNs handle well high-dimensional data and capture intricate, non-linear, and complex relationships by learning representations that fall within the graph’s topology, they are particularly, suited for detecting latent factors, thereby providing insights into the underlying structure that traditional methods might miss. This is due to the internal iterative and dynamic training and validation processes built into the GNN models. This section puts into play these strengths of GNN models to unearth latent factors of longitudinal panel data.

This section extends a model called Adversarially Regularized Graph Variational Auto-Encoder (ARGVA) architecture (Pan et al. Reference Pan, Hu, Long, Jiang, Yao and Zhang2018), to extract latent space features conditioned with the TDA features. The two key aspects of this model are for being adversarial and a variational auto-encoder, plus the important role that regularization plays overall. A generative GNN model can additionally be specified by imposing a topological constraint on extracting latent variables by adding a regularizer to the process of propagation through the various neural networks. This is based on the argument that a latent variable must be conditioned by the deep seated topological properties of the observed data since it is generated in the same manifold as the observed data, the deep properties of which are reflected in the TDA-unearthed topological factors. The model is schematically represented in Figure 8, where $\zeta $ represents the topological effects constraint and Z the latent variables.

Figure 8

Conditioned ARGVA model with the topological conditioning done through embedding: $Z \times \boldsymbol{\zeta } \to Z.$

The model effectively encodes the structure, node, and edge contents of a graph into a compact representation-embedding. The embedding is moderated by the topological effects ( $\zeta $ ). The embedding is then submitted to a decoder which is trained to reconstruct the input graph structure. The intermediate latent representation is internally forced to match a learnable probability distribution through an adversarial training (discriminator) module. We then jointly optimize the graph encoder learning and adversarial regularization to obtain the best graph embedding in a lower dimensionality compact and continuous feature space (from a latent space with, for example, 64 dimensions to a two-dimensional space in this article because the illustration has two features per node/state), while preserving the information about the graph structure, topological constraint, and node and edge features into the embeddings (Pan et al. Reference Pan, Hu, Long, Jiang, Yao and Zhang2018; Zhang et al. Reference Zhang, Yin, Zhu and Zhang2017).Footnote ¹⁰ I use EdgeGATConv as the base layer in the encoding, decoding, and discriminating phases, which includes both node and edge features in the computation as well as attention mechanisms. The Encoder, Decoder, and Discriminator qua neural networks are constituted of large numbers (100) of these layers stacked together, with thousands of learnable parameters. Formally, the neural networks are instantiated through so-called weights, which is what is learned in the training (see Supplementary Material for mathematical details) and then saved. Once the latent space is learned, I use PCA to reduce the high-dimensionality space into two principal components. To preserve the structure of the original data, the model learns separately the latent space features for each graph (year). We are working not with whole graphs (as in the classification model of the previous section) but rather with data within each graph as the task is to reconstruct each graph data separately to generate latent variables that are faithful for every panel of the data. I, thus, split the nodes and corresponding edges into training, validation, and testing subsets in each graph.

Figure 9 displays the KDE (Kernel Density Estimation) for year 2010 (as an illustration) for each for the two possible values (0,1) of the MID variable disp for training and testing data. The overall shapes for both latent space features are not that different going from training to testing data for both values (0,1) of the MIDs. The variability in the latent factors between the two values (0,1) of the MIDs suggests a possible statistical correlation between MID and the latent factors. This will be probed using a Bayesian model as done in previous sections.

Figure 9

KDE of latent space PCA components in year 2010.

As a second illustration, consider the question of nodes (states) memberships in the topological classes (clusters), and how it changes over the years. We can plot the time variation of the cohesion of all clusters, where cohesion indicates how close the nodes within the same cluster are to each other. We can then compare its trending behavior, for example, to the variation of MIDs, as shown in Figure 10, which displays spline-smoothed plots of both the cohesion measures of the clusters and the MIDs measure over time. The plot provides a visualization of the dynamic relationship between international system cohesion and military disputes. It suggests that while global disputes increased with rising cohesion during the mid-20th century, the relationship became more complex towards the end of the century, with some clusters maintaining or increasing cohesion even as overall disputes declined. During the 1946–1965 period, three clusters (0,2,3) show increasing cohesion, which coincides with an increase in the number of MIDs. Cluster 1 remains almost constant until about 1980. This suggests that as the international system became more cohesive, possibly due to alliances or blocs at the beginning of the Cold War and as many formerly colonized states achieved independence and leaned toward either side of the Cold War. It is also the era when conflicts intensified such as in the Korean war, Vietnam war, and a number of wars of independence. From 1985 to 2010, after the peak of the Cold War, the number of MIDs declines, and we observe varied cohesion trends across clusters. Three clusters (1,2,3) display a maximum and then decline between 2000 and 2010, whereas Cluster 0 reaches a local minimum and then bounces back upward. The decline in MIDs matches the behavior of clusters (1,2,3). This corresponds to the end of the Cold War, followed with a decrease in large-scale international disputes. The differing cohesion trends seem to indicate different regional dynamics or new forms of international cooperation, or conflict.

Figure 10

Cluster cohesion over time.

Figure 11 displays the cluster membership probabilities for six states around the world. The USA and RUS (Russia) plots show relatively stable membership probabilities across different clusters, indicating their central roles in global systems that persist over time. For instance, after 1960 the USA remains consistently likely to belong to all four clusters, with some variation. Similarly, Russia has high membership probabilities in Clusters 0 and 2, with a noticeable transition period in the early 1990s, likely reflecting the post-Soviet domestic and geopolitical changes. The China plot displays a more dynamic shift between clusters, especially after the early 1970s, which correspond to the period around the People’s Republic of China became more active internationally during the “rapprochement” with the US. The UK, NIG (Nigeria), and BRA (Brazil) plots exhibit more fluctuation in cluster membership, indicating that their roles or the nature of their international interactions are more variable or influenced by regional and global changes. Nigeria and Brazil, for example, have periods where they shift significantly rapidly between different clusters, likely reflecting changes in their respective regional influences and external relations.

Figure 11

Probabilities of state membership in four clusters.

Building on these results, I now consider a Bayesian model that includes all empirical variables, the expanded SAR variables, the topological effects, and the latent variables. The latter are construed as moderators much like the topological effects (with a similar specification in the model equation). The modified model specification of the augmented SAR model Equation (5.1) is as follows:

(5.1) $$ \begin{align} \gamma_i &= \sum_{j=1}^{n_{\text{vars}}} X_{ij} \left( \sum_{k=1}^{n_H} H_{ik} \delta_k \right) + \sum_{k=1}^{n_H} H_{ik} \zeta_k , \textsf{ topological moderation effects } \nonumber\\ \omega_i &= \sum_{j=1}^{n_{\text{vars}}} X_{ij} \left( \sum_{k=1}^{n_L} L_{ik} \theta_k \right) + \sum_{k=1}^{n_L} L_{ik} \xi_k , \textsf{ latent-space moderation effects } \nonumber\\ \eta_i & = \gamma_i + \omega_i + \left(1 + \rho \mathbf{W} + \rho^2 \mathbf{W^2} + \rho^3 \mathbf{W^3} + \rho^4 \mathbf{W^4} \right)\cdot \mathbf{X}_i \cdot \boldsymbol{\beta} \nonumber \\ y_i &\sim \textsf{ZeroInflatedBinomial}(\psi, N_{\text{obs}}, \text{sigmoid}(\eta_i)), \end{align} $$

where $\boldsymbol{\xi }$ are the coefficients for latent-space effects; $\boldsymbol{\theta }$ coefficients for moderation terms; $\mathbf {L} $ is the matrix of moderators/latent factors (dimension: $ N \times n_L $ ), with the rest of the quantities defined as before in Equation (4.1). A Bayesian estimation leads to the results in Table 3.

Table 3

All encompassing model.

I compare these results and all other estimations using the ELPD measure (Vehtari et al. Reference Vehtari, Gelman and Gabry2017). The all-inclusive model and the model including the sociomatrix and topological effects are very close to one another with a 0.3% difference in the elpd_loo measure. This is not surprising as we see on Table 3 that the mean for the latent variables and moderators are very small and the HDI intervals both include zero value.Footnote ¹¹ However, they both perform much better than the bare and SAR models. Continuing the task of probing how to make use of the information about topological heterogeneities, I introduce next another GNN model used to optimally reduce the dimensionality of data.

5.4. GNN-Based Dimensionality Reduction of Data

This section proposes an approach that radically reduces the dimensionality of data, while still preserving essential properties of the original data. It is based on a GNN classification model with topological factors taken as graph labels for graphs and MIDs values as node labels.

The entry point to building this approach is an off-shoot to the question of explaining the decisions underpinning the learning process in GNN models, a very active area of research (Agarwal et al. Reference Agarwal, Queen, Lakkaraju and Zitnik2023; Zhou et al. Reference Zhou2020). SubGraphX model addresses this issue by constructing an algorithm that finds subgraphs of a graph that are most influential in a model’s predictions using a Monte Carlo Tree Search (MCTS) based strategy; that is, finding subgraphs, when removed, most significantly affect the output of the GNN model. This clearly makes it a good candidate for dimensionality reduction of complex data. This idea is formalized in Equation (5.2).

(5.2) $$ \begin{align} \text{Score}(S) = \frac{1}{N} \sum_{i=1}^{N} L(y, f(G \setminus S; \theta)). \end{align} $$

S represents a candidate subgraph, $G \setminus S$ denotes the graph G excluding the subgraph S, y is the ground truth label (which in the illustration corresponds to MID values, (0,1)), f denotes the GNN with parameters $\theta $ , and L is a loss function measuring the discrepancy between the GNN output and the ground truth. The importance of each subgraph is evaluated using Shapley values as scores.Footnote ¹² Subgraphs with higher scores are considered more influential on the output of the GNN, and, therefore, good explainers of the behavior of the model performance.

Figure 12

Most important subgraphs for four years.

To illustrate, Figure 12 displays four such subgraphs obtained for different years.Footnote ¹³ The basic layer-architecture of the GNN model that I use to generate the subgraphs is similar to the one used in the previous section. The saved trained model is deployed to find the most important subgraphs for each year. I then convert these sugraphs into a longitudinal panel data the size of which will be much smaller (1118 instead of 240,273 for the original data). The generated data is then used to estimate the basic Bayesian model defined before. The results are displayed in Table 4.

Table 4

Bayesian estimates with SubGraphX-generated data.

Figure 13

ROC curve using empirical data as testing data.

To test the robustness of this model, I run an out-of-sample predictive check using the original empirical data (minus the data for those dyads that are included in the SubgraphX-generated data to avoid overfitting). I use the ROC AUC curve to display the results shown in Figure 13. The AUC value indicates a 75% chance that the model will be able to distinguish between a random positive and a random negative example. The accuracy level is 99.4%. Given that the SubgraphX-generated longitudinal panel data has just 1118 obervations as compared to 242,000 observations of testing data, the predictive power is quite impressive.Footnote ¹⁴