A network analysis of skill-specific internal migration flows in Italy☆,☆☆,☆☆☆,☆☆☆☆

The primary abstraction of the Italian internal migration network from data is a weighted and directed network. For each annual period, we created a unique network. Within this network, the links were emblematic of individual mobility, signifying the transition of people from an origin province to a destination province. The weight assigned to each link precisely quantified the volume of migrating individuals. In a formal context, we characterized this network as $G\left(V,E,W\right)$ where $\left|V\right|$ denoted the set of nodes identified as $i$, $j$ , and so forth, $\left|E\right|$ represented the set of directed edges, and $W$ symbolized the weighted adjacency matrix. This structured network encompassed $n$ nodes and $m$ edges, with the matrix element $W_{ij}$ illustrating the count of migrants moving from node $i$ to node $j$. We exclude self-loops, so the weight between a node $i$ and itself ($W_{ii}$) is always 0, denoting, in our practical application, internal migration within the same province. The total outflow from node $i$ is denoted as $W_{i.}$, and the total inflow to node $j$ is denoted as $W_{.j}$. The sum of all weights in the network is denoted as $W_{..}$. This systematic representation provided a comprehensive understanding of the migratory dynamics over distinct years, capturing the nuances in the flow of individuals between origin and destination provinces.

In network analysis literature, there is an active field focused on extracting meaningful structures from dense networks, eliminating noise-like connections. This pursuit, known as backbone extraction, aims to pinpoint significant connections for each object within the network [34]. The implementation of this approach varies across dimensions, categorized into statistical and structural methods, with hybrid approaches combining both [35]. Statistical methods employ hypothesis testing and p-values, while structural methods work on the network topology to derive a backbone with specific characteristics.

This paper addresses backbone extraction in migration networks. Migration networks are inherently dense and characterized by heterogeneous node strengths and edge weights, indicating a multiscale nature. The goal of deriving a meaningful signed network is achieved in this work by employing the method proposed by [34]. We denote the signed network as $\hat{G}=\left(\hat{V},\hat{A},\hat{W}\right)$, where $\hat{A}$ is a sparse adjacency matrix with values $\left\{-1,0,1\right\}$, indicating negative, absent, or positive links. The corresponding weight matrix is denoted as $\widehat{W.}$ For extracting the signed edges in a network is crucial to establish a null model for comparing observed edge weights with expected values, allowing differentiation between chance occurrences and significant deviations in either a positive or negative direction. In the null model, we assume that each node in-strengths and out-strengths are predetermined, meaning that each node possesses a specific number of stubs or half edges. A unit-link signifies a connection between two stubs from different nodes. When a stub seeks connection, it randomly selects another available stub. The likelihood of connecting to a specific node is proportional to the number of available stubs for that node. The weight of the link between two nodes reflects the number of unit-links connecting them. This process resembles a sampling without replacement problem, akin to the well-known urn problem. To better understand, consider node $i$: there are $W_{..}-W_{.i}$ marbles in the urn (excluding loops), $W_{i.}$ marbles are drawn from the urn without replacement, and there are $W_{.j}$ marbles in the urn for each $j$. This process aligns with the hypergeometric distribution, facilitating the computation of statistical quantities for the links between all nodes $i$ and $j$. The mean of the hypergeometric distribution, associated with the link from node $i$ to node $j$, is denoted as $P_{ij}$ and is determined by Eq. (1). However, this formulation does not preserve the sequences of in-strength and out-strength in the system, making it inappropriate for the intended purposes.(1)$P_{ij}=\frac{\left(W_{i.}{W}_{.j}\right)}{\left(W_{..}{-W}_{.i}\right)}$(2)$\underset{N_{ij}}{\min }\sum _{i,j\ne i}{N}_{ij}\log \frac{N_{ij}}{P_{ij}}$$subjectto{N}_{i.}=W_{i.},N_{.j}=W_{.j},N_{ij}\ge 0,N_{ii}=0,\forall i,j\text{.}$

To ensure the preservation of in-strength and out-strength sequences, certain conditions must be met, as detailed in Eq. (2), where $N$ is the weight matrix under the null model and $N_{ij}$ represents the expected weight from node $i$ to node $j$ in this null model.

The iterative proportional fitting procedure (IPFP), described by Ref. [36,37], is adopted to meet the above requirements. IPFP utilizes the out-strength and in-strength sequences ($W_{i.}$ and $W_{.j}$) as marginals, and the prior matrix $P$, where its diagonal elements are set to 0. The objective of IPFP is to estimate the weight matrix $N$ that minimizes the relative entropy, also known as the Kullback-Leibler divergence [38], as formulated in Eq. (2).

IPFP iteratively adjusts the values of non-zero elements in the matrix through row-scaling and column-scaling operations until a desired level of precision is attained. It has been demonstrated that IPFP converges to an optimal solution [39]. Each element of the estimated weight matrix $N$ represents the expected value of the corresponding link under the null model being utilized. To further analyze the estimated values and their reliability, it is necessary to establish confidence intervals. These intervals rely on a dispersion measure, with the variance of $N_{ij}$ estimated using the associated hypergeometric distribution. The well-known standard deviation is then obtained from this variance, as given by Eq. (3):(3)${\sigma }_{ij}=\sqrt{\left(W_{i.}\frac{W_{.j}}{W_{..}-W_{.i}}\frac{\left(W_{..}{-W}_{.i}-W_{.j}\right)}{W_{..}{-W}_{.i}}\frac{W_{..}-W_{.i}-W_{i.}}{W_{..}-W_{.i}-1}\right)}$

A combination of two filters, the significance filter, and the vigor filter, is then applied to eliminate links with weights that do not significantly deviate from expected values. The significance filter removes links that meet the criteria,${\sigma }_{ij}{\alpha }^{-}<W_{ij}-N_{ij}<{\sigma }_{ij}{\alpha }^{+}$where ${\alpha }^{+}$ and ${\alpha }^{-}$ are user-defined hyperparameters specifying the desired significance thresholds for negative and positive signed edges. Since the significance filter might be too permissive, retaining very weak links, we also consider a vigor filter, expressed as${\beta }_{ij}=\frac{W_{ij}/N_{ij}-1}{W_{ij}/N_{ij}+1}$${\beta }_{ij}$ takes values in the range $\left[-1,1\right]$, with ${\beta }_{ij}=0$ when $W_{ij}=N_{ij}$. Vigor filter remove links that fall within the range defined by the condition that${\beta }^{-}<{\beta }_{ij}<{\beta }^{+}\text{,}$where ${\beta }^{+}$ and ${\beta }^{-}$ are user-defined hyperparameters determining the desired thresholds for positive and negative signed edges.

The resulting signed backbone, extracted using both filters, displays features commonly linked to signed networks, including reciprocity, structural balance, and community structure. It is worth noting that, extracting negative links from networks with positive edge weights relies on the assumption that links with small edge weights or the absence of links may represent negative connections. While this assumption may not apply to many networks, it is deemed acceptable in intrinsically dense networks, such as migration networks.

Efficient analysis of network data relies on the chosen representation [40]. Traditional forms, like adjacency matrices, only capture node relationships with neighbors, presenting computational challenges for large-scale networks. To overcome this, high-dimensional data is reduced to a more manageable lower-dimensional structure. Network representation learning (RL), a dimensionality reduction technique [41], achieves this by mapping high-dimensional network data to a lower-dimensional vector space while preserving network topology. It minimizes redundancy and noise, retaining essential structural details. The acquired vertex representations must meet three conditions:

• 1.
  they should be of low dimensionality, much smaller than the size of the vertex set;
• 2.
  they should be informative, maintaining the proximity of vertices through consideration of network structure, attributes, and labels;
• 3.
  they should be continuous, with real values, facilitating tasks such as vertex classification and clustering.

The essential goal is to construct an embedding matrix $Z$ with dimensions $\left(n\times d\right)$, where $n$ represents the number of nodes, $d$ is the embedding size and $z_i$ denotes the embedding vector for node $i$. This ensures that the similarity within the embedding space closely mirrors the similarity observed in the network. Formally, the task involves estimating the embedding matrix $Z$ while satisfying the condition $f\left(i,j\right)\sim g\left(z_i,z_j\right)$, where $f$ is a similarity function in the observed network, and $g$ is a similarity function in the latent space.

In our setting, we begin by determining the underlying signed backbone network $\hat{G}=\left(\hat{V},\hat{W}\right)$, with $\hat{W}$ representing the similarity in the observed network. For measuring similarity in the latent space, we opt for cosine similarity, characterized by a natural range of $\left[-1,1\right]$, aligning with the range of the input data. Once we define the similarity functions in both the observed space and latent space, we proceed to formulate the loss function $\mathcal{L}$ as outlined in Equation (4).(4)$\sum \left|cossim\left(i,j\right)-W_{ij}\right|$

Stochastic gradient descent [42] is employed to discover the optimal $Z$ that minimizes $\mathcal{L}$. The process involves initializing $Z$ with random values, iterating over non-diagonal elements of $W$ in a randomized order, calculating gradients for $z_i$ and $z_j$, and updating them using a learning rate $\eta =0.01$. Each iteration over all elements of $W$ is termed an epoch, and this process is repeated for 100 epochs for each network in a given year.

The optimization process may converge to diverse yet equally plausible solutions, stemming from the random initialization of $Z$ or the intrinsic stochastic nature of the optimization procedure. Notably, any rotation or reflection of the determined solution yields identical information. Consequently, a direct comparison of embedding vectors between different time steps is not feasible. Optimal rotations and/or reflections might exist between embeddings of different years, preserving cosine similarity while altering the spatial orientation. This occurrence is identified as the Procrustes problem [43], extensively discussed by Gursoy & Badur [34].

In the last step, we address the objective of partitioning network vertices into clusters, where vertices within the same cluster display dense connections while having fewer connections to vertices from other clusters [44]. These cluster structures, often referred to as communities, are prevalent in a diverse range of networked systems, and carry significant implications. In migration networks, clusters represent provinces with concentrated migration flows. Analyzing spatial interactions aids in understanding dependencies, cultural exchanges, and migration's impact on interconnected areas, guiding policy decisions and resource allocation while predicting future migration trends. The literature offers a variety of network clustering algorithms that utilize different metrics to measure similarity or the strength of connections between vertices [[45], [46], [47], [48]]. In this context, to understand the community structures in the network, we employed the agglomerative hierarchical clustering [49], based on learned representations to cluster vertices and density-based spatial clustering of applications with noise (DBSCAN [50]). DBSCAN groups data points based on density, forming clusters and marking outliers. Its effectiveness relies on hyperparameters, such as eps determining neighborhood radius and minSamples specifying the minimum objects for a dense region. Adjusting these parameters is crucial, as smaller eps values create denser clusters, while larger values may merge or introduce noise. Conversely, elevating minSamples values generally yields larger and more stable clusters, whereas smaller values can result in more clusters and heightened sensitivity to noise.

In the hierarchical method, we measure the dissimilarity between two clusters by considering the maximum distance between their respective data points, using the complete linkage method [49].

The distance calculation between two provinces, $i$ and $j$, is expressed as $1-cossim\left(z_i,z_j\right)$, where cosine similarity quantifies the similarity in the latent space. Consequently, we derive distances for all province pairs within the $\left[0,2\right]$ range. These distance values serve as inputs for our clustering analysis.