Generalized gravity model for human migration

For decades, the gravity model (GM) has successfully explained flows between geographically separated two regions such as traffic flow [1–5], international economic trades [6, 7], and human migration [8, 9]. The GM is named after Newton's law of universal gravitation because of the similarity in the formula: a certain type of flow between two regions is proportional to the product of 'mass' of each region and inversely proportional to a certain power of distance between the regions. We interpret the mass depending on contexts; we can quantify the relative importance of regions in human migration from their population sizes, and the relative importance of countries in international trades from their economic scales. This simple but powerful model has succeeded in interpreting the real world. For instance, the GM indeed accurately describes the empirical data of daily human mobility in multiscale mobility networks [10]. It also nicely explains the inter- and intra-city traffic flows in Korea [11], along with the passenger flows in the Korean subway system [12].

However, those examples only concern the spatial aspect of population. We can easily imagine more complicated situations such as the population flow between other attributes than the spatial or regional attributes, when the population at a given region consists of subpopulation structures. The subpopulation structures characterized by attributes for population flows can be ethnic groups, income levels, etc. In particular, when spatial movement of the subpopulation belonging to one attribute to another subpopulation takes place, the flow between attributes becomes relevant. This applies not only to the population flow but also to the international trade, for instance. Goods are transferred in different economic sectors as the attributes, and one may aim at estimating the flow between economic sectors. One such previous attempt to apply the GM to explain flow between attributes is [13] partly by a subset of the authors of this paper. However, the results have revealed the limitation of the conventional GM when the geographical location of population center of a clan is not representative. The limitation stems from the process of significant coarse graining of the detailed geographical information of clans into a single point (the population center).

In this paper, we exploit detailed information on the population substructures by generalizing the GM, rather than coarse graining the information. First, we formulate the generalized gravity model (GGM) to take the full information available on the subpopulation (characterized by both geographical information and attribute information) flows. As we show in figure 1(a), when subpopulations are distinguished not only by geographical regions but also by attributes, flows between regions or attributes can be calculated by properly taking subpopulation flows without information loss (see figures 1(b) and (c), respectively). We also would like to emphasize the necessity of calculating subpopulation flows, because the population flows calculated from the coarse-grained population data and those from individual subpopulation data are not equivalent. As a concrete example, we apply the GGM to marriage records combined with census data. The results show that it effectively captures the geographical constraint imposed in the marriage patterns in the past, in contrast to the GM.

Zoom In Zoom Out Reset image size

The paper is organized as follows. We first formulate the GGM in section 2. We introduce our data set in section 3 and apply the model to this data. The result in section 4 demonstrates that the GGM indeed captures the geographical information not available from the GM. Finally, we conclude our work in section 5.

Our derivation of the GGM is a natural extension of using the maximum entropy principle to derive the GM [14, 15], where we replace regional indices with both regional and attribute indices. We provide our step-by-step derivation in this section partly for a pedagogical purpose and the self-containedness of this paper, but most importantly, we can directly demonstrate the problem of using the coarse-grained population data during the derivation. The maximum entropy principle is the way to estimate a real probability distribution by maximizing entropy. In particular, this method is useful for systems with many degrees of freedom because it focuses on only a few macroscopic quantities. The real probability distribution is estimated by the maximum entropy principle from the agreement of those observed quantities. Each Lagrange multiplier corresponding to each constraint in maximization gives the corresponding model parameter.

Let us start from the number ${N}_{{AB}}^{({ij})}$ of people who move from region A with attribute i to region B with attribute j, where we take the convention of uppercase letters as the subscript for the region indices, and the lowercase letters as the superscript for the attribute indices. The sets of attributes {i} for the sender side and {j} for the receiver side can be different. For example, {i} and {j} can be the education level and the income level, respectively, when we try to describe the flow of employment from one city's education system to another city's industry. Our formalism adopts the discrete indices and summation, but it should be straightforward to deal with the continuous cases by using continuous variables and integration.

The total number ${N}_{A}^{(i)}$ of people moving from region A with attribute i to anywhere with any attribute is then

$$\begin{eqnarray}&&{N}_{A}^{(i)}=\displaystyle \sum _{j,B}{N}_{{AB}}^{({ij})}.\end{eqnarray} \tag{ 1 }$$

In the same way, the total number ${\tilde{N}}_{B}^{(j)}$ of people who arrive at region B with attribute j from anywhere with any attribute is given by

$$\begin{eqnarray}&&{\tilde{N}}_{B}^{(j)}=\displaystyle \sum _{i,A}{N}_{{AB}}^{({ij})}.\end{eqnarray} \tag{ 2 }$$

When people move, they have to pay the cost, which is naturally a function of the distance between two regions, among other factors. We denote the cost to move from (i, A) to (j, B) for each unit of movement by ${c}_{{AB}}^{({ij})}$. Then, the total moving cost C is the following weighted sum

$$\begin{eqnarray}&&C=\displaystyle \sum _{i,A,j,B}{N}_{{AB}}^{({ij})}{c}_{{AB}}^{({ij})}.\end{eqnarray} \tag{ 3 }$$

We can then write down the number W of all possible arrangements of travelers considering the multiplicity factor ${N}_{{AB}}^{({ij})}$ as

$$\begin{eqnarray}&&W=\displaystyle \frac{N!}{{\prod }_{i,A,j,B}{N}_{{AB}}^{({ij})}},\end{eqnarray} \tag{ 4 }$$

where N is the total number of moving people, $N={\sum }_{i,A,j,B}{N}_{{AB}}^{({ij})}$. In the entropy maximization scheme, $\{{N}_{{AB}}^{({ij})}\}$ is estimated from maximizing the Boltzmann entropy ${k}_{B}\mathrm{log}W$ (equivalently maximizing W) with constraints. We consider three constraints under which W is maximized: the outflows $\{{N}_{A}^{(i)}\}$, the inflows $\{{\tilde{N}}_{B}^{(j)}\}$, and the total moving cost C represented in equations (1)–(3). For this optimization problem under given constraints, we have to use the Lagrange multiplier method, i.e., to find the stationary point of the Lagrangian

$$\begin{eqnarray}&&{ \mathcal L }(\{{N}_{{AB}}^{({ij})}\})=\displaystyle \frac{N!}{{\displaystyle \prod }_{i,A,j,B}{N}_{{AB}}^{({ij})}}+\displaystyle \sum _{i,A}{\lambda }_{A}^{i}\left[{N}_{A}^{(i)}-\displaystyle \sum _{j,B}{N}_{{AB}}^{({ij})}\right]+\displaystyle \sum _{j,B}{\tilde{\lambda }}_{B}^{j}\left[{\tilde{N}}_{B}^{(j)}-\displaystyle \sum _{i,A}{N}_{{AB}}^{({ij})}\right]+\gamma \left[C-\displaystyle \sum _{i,A,j,B}{N}_{{AB}}^{({ij})}{c}_{{AB}}^{({ij})}\right],\end{eqnarray} \tag{ 5 }$$

where ${\lambda }_{A}^{i}$, ${\tilde{\lambda }}_{B}^{j}$, and γ are the Lagrange multipliers for each constraint. This problem is essentially the recap of the derivation for the most probable distribution in terms of a given energy value, from the standard formalism of canonical ensemble in statistical mechanics, so the readers may check the details in any standard statistical mechanics textbooks such as [19]. The moving cost for each movement and the γ parameter here play the roles of energy and inverse temperature there, respectively. The solution of maximizing equation (5) is given by

$$\begin{eqnarray}&&{N}_{{AB}}^{({ij})}\propto {N}_{A}^{(i)}{\tilde{N}}_{B}^{(j)}{{\rm{e}}}^{-\gamma {c}_{{AB}}^{({ij})}}.\end{eqnarray} \tag{ 6 }$$

Note that all ${\lambda }_{A}^{i}$ and ${\tilde{\lambda }}_{B}^{j}$ become unity from the constraint itself (the mass conservation), so there is only one free parameter γ, which is determined by the real data. Later, we will specifically choose the γ value that minimizes the error between the model and real data.

The flow from region A to region B usually decays as a function of the distance between them due to the obviously rising cost, and thus we have to choose the cost function c(r) as an increasing function of the distance r between two regions. Conventionally, we set the form $c(r)\propto \mathrm{ln}r$ [20], which leads equation (6) to

$$\begin{eqnarray}&&{N}_{{AB}}^{({ij})}\propto \displaystyle \frac{{N}_{A}^{(i)}{\tilde{N}}_{B}^{(j)}}{{({r}_{{AB}})}^{\gamma }}.\end{eqnarray} \tag{ 7 }$$

Note that, though the numbers of leaving or arriving people, ${N}_{A}^{(i)}$ and ${\tilde{N}}_{B}^{(j)}$, are not the population size ${P}_{\mathrm{region}}^{\mathrm{attribute}}$ at each region and attribute, generally those are assumed to be linear to the population sizes [${N}_{A}^{(i)}\propto {P}_{A}^{(i)}$ and ${\tilde{N}}_{B}^{(j)}\propto {P}_{B}^{(j)}$], yielding

$$\begin{eqnarray}&&{N}_{{AB}}^{({ij})}\propto \displaystyle \frac{{P}_{A}^{(i)}{P}_{B}^{(j)}}{{({r}_{{AB}})}^{\gamma }}.\end{eqnarray} \tag{ 8 }$$

In this case, we can reproduce the GM for the flow between two spatial regions

$$\begin{eqnarray}&&{N}_{{AB}}\equiv \displaystyle \sum _{i,j}{N}_{{AB}}^{({ij})}\propto \displaystyle \frac{{\sum }_{i}{P}_{A}^{(i)}{\sum }_{j}{P}_{B}^{(j)}}{{({r}_{{AB}})}^{\gamma }}=\displaystyle \frac{{P}_{A}{P}_{B}}{{({r}_{{AB}})}^{\gamma }},\end{eqnarray} \tag{ 9 }$$

where P_A and P_B are total numbers of people who live in region A and region B, respectively. However, when ${N}_{A}^{(i)}$ and ${\tilde{N}}_{B}^{(j)}$ are nonlinear with respect to the population sizes such as in [10], i.e., ${N}_{A}^{(i)}\propto {[{P}_{A}^{(i)}]}^{\alpha }$ and ${\tilde{N}}_{B}^{(j)}\propto {[{P}_{B}^{(j)}]}^{\beta }$, the population flow from region A with attribute i to region B with attribute j: ${N}_{{AB}}^{({ij})}\propto {[{P}_{A}^{(i)}]}^{\alpha }{[{P}_{B}^{(j)}]}^{\beta }/{({r}_{{AB}})}^{\gamma }$. In this case, the flow from region A to region B regardless of attributes: N_AB cannot be explained by the GM unless α = β = 1 or the population contains only one attribute, because

$$\begin{eqnarray}&&{N}_{{AB}}\equiv \displaystyle \sum _{i,j}{N}_{{AB}}^{({ij})}\propto \displaystyle \frac{{\sum }_{i}{[{P}_{A}^{(i)}]}^{\alpha }{\sum }_{j}{[{P}_{B}^{(j)}]}^{\beta }}{{({r}_{{AB}})}^{\gamma }}\ne \displaystyle \frac{{[{P}_{A}]}^{\alpha }{[{P}_{B}]}^{\beta }}{{({r}_{{AB}})}^{\gamma }}.\end{eqnarray} \tag{ 10 }$$

Note that the conventional GM provides the upper or lower bounds: ${\sum }_{i}{[{P}_{A}^{(i)}]}^{\alpha }\leqslant {[P]}^{\alpha }$ for α > 1 and ${\sum }_{i}{[{P}_{A}^{(i)}]}^{\alpha }\geqslant {[P]}^{\alpha }$ for α < 1, by using the convexity or concavity of the functional form.

In parallel, summing up ${N}_{{AB}}^{({ij})}$ for all of the regions gives the number of moving people from attribute i to j,

$$\begin{eqnarray}&&{N}^{({ij})}\equiv \displaystyle \sum _{A,B}{N}_{{AB}}^{({ij})}\propto \displaystyle \sum _{A,B}\displaystyle \frac{{N}_{A}^{(i)}{\tilde{N}}_{B}^{(j)}}{{({r}_{{AB}})}^{\gamma }}.\end{eqnarray} \tag{ 11 }$$

This is not reducible to the GM either, because

$$\begin{eqnarray}&&\displaystyle \sum _{A,B}\displaystyle \frac{{N}_{A}^{(i)}{\tilde{N}}_{B}^{(j)}}{{({r}_{{AB}})}^{\gamma }}\ne \displaystyle \frac{1}{{[r({i}_{\mathrm{CM}},{j}_{\mathrm{CM}})]}^{\gamma }}\displaystyle \sum _{A,B}{N}_{A}^{(i)}{\tilde{N}}_{B}^{(j)},\end{eqnarray} \tag{ 12 }$$

where r(i_CM, j_CM) is the distance between the centers of population of i and j, as shown in figure 1. The only case when the two expressions actually coincide is the assumption implicitly made in [13]—the population of each attribute i is treated as the 'point mass' located in a single location in space, namely, i_CM. The GGM estimates the flow between attributes without such a coarse graining process involving the information loss. Hence, the GGM is the correct way to handle subpopulation structures when it comes to the GM of population flow. We later show that it indeed effectively captures the geographical constraints for the marriage flow in the past obtained from the data, which was not possible with the coarse-grained version of the GM due to the widely distributed population of the clans [13].

In the traditionally patriarchal culture of Korea after around 17th century, a bride usually moved to her groom's place in the past, once they got married. We treat this type of migration caused by the marriage as our main data of human migration. By applying the GGM to marriage patterns between clans in the past, we estimate the geographical constraints. We take the real marriage flow ${O}^{({ij})}$ from the bride clan i to the groom clan j from the family book data called jokbo. We present more details in section 3.1. To compute the model flow, we extract the distance r_AB between two regions and the distribution of the population for each clan, ${N}_{A}^{(i)}$ (the bride side) and ${\tilde{N}}_{B}^{(j)}$ (the groom side), from the modern census data in 1985, 2000, and 2015—the three particular years when the information on the regional distribution of each clan's population is available.

We measure the distance r_AB based on geographical coordinates of the regions using the Google maps application programming interface [21]. The traveling distance within the same region r_AA is estimated as the square root of the region's area. We assume that the size of the moving populations, represented by ${N}_{A}^{(i)}$ and ${\tilde{N}}_{B}^{(j)}$, are proportional to that of the resident populations (from the census data) of the corresponding clans living in the corresponding regions, so we just take the face values of populations in the census data and treat them as the migrating population for simplicity. As argued in [13], we use this modern population data to estimate the past migration flow between clans in jokbo data, based on the fact that the proportion of each clan living in each region with respect to the total population of Korea has been relatively steady.

3.1. Jokbo data

3.2. Census data

We apply the GGM expressed in equation (11) to the marriage patterns in the past. The real number of marriage entries {${O}^{({ij})}$} are listed in the jokbo data, and we compare the predicted flow $\{{N}^{({ij})}\}$ from the model with {${O}^{({ij})}$}. For each groom clan (corresponding to each jokbo clan) j, the difference is quantified by the error

$$\begin{eqnarray}&&{E}^{(j)}=\sqrt{\displaystyle \sum _{i}{[{E}^{({ij})}]}^{2}}\equiv \sqrt{\displaystyle \sum _{i}{[{O}^{({ij})}-{N}^{({ij})}]}^{2}},\end{eqnarray} \tag{ 13 }$$

calculated from the list of the bride clans {i}. Note that we discard the self migration flow ${E}^{({jj})}$ when we calculate ${E}^{(j)}$, because the marriage between the same clans was forbidden in the past and it is indeed significantly underrepresented as reported in [13]. In practice, we also checked that ignoring ${E}^{({jj})}$ does not make much of a difference in our results. The proportionality factor for equation (11) is calculated by minimizing ${E}^{(j)}$ at a given value of ${\gamma }^{(j)}$. The optimal value ${\gamma }_{\mathrm{opt}}^{(j)}$ is assigned as the ${\gamma }^{(j)}$ value that minimizes the error ${E}^{(j)}$. In this case, we vary the ${\gamma }^{(j)}$ value from 0 to 10 with the resolution of 0.1. The obtained ${\gamma }_{\mathrm{opt}}^{(j)}$ value indicates the geographical constraint for the brides' migration to the groom clan j.

In figure 2, on the left side of each panel (corresponding to each groom clan j), we show the error in equation (13) as a function of ${\gamma }^{(j)}$, via the fact that ${N}^{({ij})}$ is a function of ${\gamma }^{(j)}$. On the right side of each panel, we also show scatter plots comparing the real flow ${O}^{({ij})}$ versus the predicted flow ${N}^{({ij})}$ from the GGM at a given ${\gamma }_{\mathrm{opt}}^{(j)}$ value with the guideline corresponding to ${O}^{({ij})}={N}^{({ij})}$. The predicted flow ${N}^{({ij})}$ and the number of entries ${O}^{({ij})}$ in jokbo are indeed close to the ${O}^{({ij})}={N}^{({ij})}$ line. Since the results from 1985, 2000, and 2015 census data are qualitatively the same, we only show the results in 2015. Except for the clan 1, the GGM actually captures nonzero ${\gamma }_{\mathrm{opt}}^{(j)}$ values, while the exponent ${\gamma }_{\mathrm{opt}}^{(j)}$ always vanishes when we use the GM for all of the jokbo [13].

Zoom In Zoom Out Reset image size

We calculate ${N}^{({ij})}$ from each census data in 1985, 2000, and 2015, so we obtain the three ${\gamma }_{\mathrm{opt}}^{(j)}$ values corresponding to each year, for each groom clan (see figure 3). From the similarity of error landscapes in each census data, the ${\gamma }_{\mathrm{opt}}^{(j)}$ values for different years are not much different for a given groom clan j (check table 2 for details), which indicates that the results of ${\gamma }_{\mathrm{opt}}^{(j)}$ are temporally robust for a given groom clan. It is also interesting to note that many ${\gamma }_{\mathrm{opt}}^{(j)}$ values are around 2 corresponding to the same formula with demographic gravitation introduced in [2]. Most importantly, compared with the GM results, i.e., ${\gamma }_{\mathrm{opt}}^{(j)}\simeq 0$ for all of the groom clans [13], the GGM indeed yields nonzero ${\gamma }_{\mathrm{opt}}^{(j)}$ values except for the clan 1. It implies that the GGM actually captures the information of the geographical constraint on the flow, in contrast to the GM where it is hard to capture this geographical information due to the coarse graining, i.e., treating the population of a clan as a point particle at a center of mass. On the contrary, the GGM uses much more detailed information of the subpopulation structure that eventually leads us to capture the actual geographical constraint imposed on the flow.

Zoom In Zoom Out Reset image size

Table 2.  The results of ${\gamma }_{\mathrm{opt}}^{(j)}$ for each groom clan j, from the 1985, 2000, and 2015 census data, respectively. For the comparison between the GGM and the GM, we provide the values of ${\rm{\Delta }}{e}^{(j)}$ (%) in equation (14) as the percentage, representing the relative performance of the GGM. We also characterize the population distribution of each clan by measuring the dispersion ${\rm{\Delta }}{R}^{(j)}$ (km) in equation (15) and the effective number of occupied regions ${n}^{(j)}$ in equation (16).

Groom clan j	Year 1985				Year 2000				Year 2015
	${\gamma }_{\mathrm{opt}}^{(j)}$	${\rm{\Delta }}{e}^{(j)}$	${\rm{\Delta }}{R}^{(j)}$	${n}^{(j)}$	${\gamma }_{\mathrm{opt}}^{(j)}$	${\rm{\Delta }}{e}^{(j)}$	${\rm{\Delta }}{R}^{(j)}$	${n}^{(j)}$	${\gamma }_{\mathrm{opt}}^{(j)}$	${\rm{\Delta }}{e}^{(j)}$	${\rm{\Delta }}{R}^{(j)}$	${n}^{(j)}$
1	0.1	0.02	151.35	104.23	0.2	0.17	148.54	81.43	0.2	0.38	146.32	76.65
2	1.7	9.98	118.39	43.15	1.7	8.50	118.43	43.18	1.7	7.57	117.83	47.21
3	1.6	17.81	149.91	71.31	1.7	18.99	142.11	62.44	1.8	16.81	139.23	60.14
4	1.7	15.01	131.46	85.29	1.7	15.17	125.76	67.34	1.8	13.51	125.52	64.48
5	1.3	0.52	143.01	54.98	1.8	1.28	146.06	51.84	2.0	2.02	146.76	62.44
6	1.7	9.43	141.09	60.74	2.0	10.12	137.44	50.50	2.1	9.24	134.61	51.60
7	5.3	6.05	155.77	91.37	5.6	9.69	152.98	74.21	6.2	10.32	149.64	71.26
8	6.0	2.88	140.40	63.56	6.6	6.64	136.73	61.83	6.9	7.08	133.10	57.67
9	7.2	2.18	139.94	87.32	7.3	2.54	136.68	70.49	7.7	2.40	134.97	69.69

To validate our model, we measure the performance of our model compared with the GM using the normalized reduced error with respect to the case of no geographical constraint, i.e., γ = 0, defined as

$$\begin{eqnarray}&&{\rm{\Delta }}{e}^{(j)}=\displaystyle \frac{[{E}^{(j)}(\gamma =0)-{E}^{(j)}(\gamma ={\gamma }_{\mathrm{opt}}^{(j)})]}{{E}^{(j)}(\gamma =0)}.\end{eqnarray} \tag{ 14 }$$

The normalized reduced error ${\rm{\Delta }}{e}^{(j)}$ quantifies the improvement of performance by using the GGM compared with the GM, which results in ${\gamma }_{\mathrm{opt}}^{(j)}\simeq 0$ for all of the clans [13]. Large values of ${\rm{\Delta }}{e}^{(j)}$ indicate significance of geographical constraints in the migration flow. Except for clans such as j = 1, 5, and 9, the normalized reduced error ${\rm{\Delta }}{e}^{(j)}$ ≃ 10%, as shown in figure 4. To demonstrate the statistical significance of geographical information of the data, we shuffle regional indices for the groom clan j to obtain the corresponding surrogate ${\gamma }_{\mathrm{opt}}^{(j)}$ values, also shown in figure 4.

Zoom In Zoom Out Reset image size

For the exceptional cases of j = 1, 5, and 9, we suspect the lack of geographical information in the data itself, as we argue. To test the statistical significance of spatial correlation in the data, we shuffle the regional indices to scramble geographical information. We examine the result of our model in this shuffled data, which is shown in figure 4. It supports that the GGM extracts more geographical information than the GM, by capturing the nonzero γ exponent. If the shuffled data gives similar results to those from original data, the original data contains a small amount of geographical information. As shown in figure 4, this situation precisely happens for the clans j = 1, 5, and 9, whose original and shuffled data give similar results. In other words, the data corresponding to those three clans originally contain less geographical information than the other clans. The small ${\rm{\Delta }}{e}^{(j)}$, therefore, does not come from the GGM but from the clan data itself. Hence, we conclude that as long as the data has enough geographical information, the GGM effectively extracts the corresponding information.

As mentioned above, geographical information of the population distribution is closely related to the model performance and the statistical significance of the nonzero ${\gamma }_{\mathrm{opt}}^{(j)}$ values. To quantify the geographical information in the distribution of clans more systematically, we introduce two measures: the dispersion that quantifies how strongly localized the populations are, and the homogeneity that focuses on how uniformly the populations occupy distinct regions. For the latter, we use the concept of the effective number of occupied regions based on the Rényi entropy for a given probability distribution, as in [23]. We define the dispersion ${\rm{\Delta }}{R}^{(j)}$ from the centroid ${{\bf{R}}}^{(j)}$ of the clan j by taking population fractions as weights:

$$\begin{eqnarray}&&{\rm{\Delta }}{R}^{(j)}=\sqrt{\displaystyle \sum _{A}{f}_{A}^{(j)}{| | {{\bf{r}}}_{A}-{{\bf{R}}}^{(j)}| | }^{2}},\end{eqnarray} \tag{ 15 }$$

where r_A is the location of the administrative region A. The population fraction ${f}_{A}^{(j)}$ is the population of clan j living in A divided by the total population of clan j, and $| | \cdots | |$ is the Euclidean norm. The population centroid of clan j is then ${{\bf{R}}}^{(j)}={\sum }_{A}{f}_{A}^{(j)}{{\bf{r}}}_{A}$. As the concept of moment of inertia or radius of gyration [13], ${\rm{\Delta }}{R}^{(j)}$ measures how (geographically) widely a certain clan is distributed. The effective number of occupied regions is defined as the reciprocal of the heterogeneity quantified by the second moment of the population fraction, given by

$$\begin{eqnarray}&&{n}^{(j)}=\displaystyle \frac{1}{{\sum }_{A}{[{f}_{A}^{(j)}]}^{2}}.\end{eqnarray} \tag{ 16 }$$

In contrast to ${\rm{\Delta }}{R}^{(j)}$, ${n}^{(j)}$ measures how many of distinct regions (regardless of their geographical location) a certain clan occupies effectively; there are scaling relations for extreme cases: ${n}^{(j)}$≃ the total number of administrative regions when the clan j is uniformly distributed to the entire set of administrative regions, while ${n}^{(j)}$ ≃ 1 when the clan j is almost exclusively living in a single particular administrative region [23].

The dispersion ${\rm{\Delta }}{R}^{(j)}$ and the number ${n}^{(j)}$ of occupied regions are usually positively correlated. However, ${\rm{\Delta }}{R}^{(j)}$ can be large even when ${n}^{(j)}$ is small, e.g., when the clan has multiple localized residential regions. Hence, we use both measures for more accurate identification of the population distributions. For all of the 788 clans in the 2015 census data, we measure both ${\rm{\Delta }}{R}^{(j)}$ and ${n}^{(j)}$ and present them as the density plot in figure 5. Among all of the clans listed in the census data, all of the groom clans corresponding to the jokbo data have relevantly large values of dispersion and the effective number of occupied regions. Note that the combination of large ${\rm{\Delta }}{R}^{(j)}$ and small ${n}^{(j)}$ (as discussed before) is observed indeed, while the combination of small ${\rm{\Delta }}{R}^{(j)}$ and large ${n}^{(j)}$ does not appear, as shown in figure 5. This contrast hints the existence of such multiple localized residential regions, which is also discussed in [13].

Zoom In Zoom Out Reset image size

Finally, we compare the relative performance of GGM with these two measures, presented in figure 6. There is a trend of increasing ${\rm{\Delta }}{e}^{(j)}$ when the population distribution has more geographical information, low dispersion and low effective number of occupation regions. The Pearson correlation coefficients between ${\rm{\Delta }}{e}^{(j)}$ and the dispersion ${R}^{(j)}$ or the effective number ${n}^{(j)}$ of occupied regions are −0.270 and −0.213, respectively, implying the anticorrelation. This observation again confirms that when the data includes meaningful geographical information, the GGM captures the geographical constraint, while the GM may not.

Zoom In Zoom Out Reset image size

In this paper, we formulate the GGM to properly take subpopulation structures for human migration, by keeping the entire geographical distribution of the subpopulations. The key aspect of our point is that we need to calculate individual subpopulation flows, before trying any geographical coarsening. To test the validity of the GGM, we investigate the marriage patterns of Korea in the past. Applying our model to the marriage pattern, we identify the geographical constraint. The results demonstrate that the GGM captures the subpopulation aspect of the data without the information loss occurred in the GM.

We believe that our approach is applicable to a wide range of research on population dynamics. Moreover, we would like to point out that the GGM is in fact even more general than the treatment for our particular data set, e.g., the different types of attributes for the departure and arrival places by taking the different sets $\{i\}\ne \{j\}$. For instance, the attribute in the departure place for education can be the education level of people, while the attribute of the arrival place for work can be the income level. Furthermore, if we release the constraints α = 1 and β = 1, one can also allow the nonlinear mass relation. In this case, the GGM becomes particularly important to prevent the loss of information because the scale factor, in addition to the distance factor, also has the nonlinear relation with the flow.

Finally, this scheme can be extended for multiple types of attributes that can be represented by the attribute vectors. For instance, people living in the region A with the attribute i = (e_i, m_i) representing the education level e_i and the income level m_i can move to the region B with j = (e_j, m_j). We hope to extend this type of general scheme for a wide variety of different data sets of human (and possibly nonhuman) migration or flow patterns in the future.

SHL was supported by Gyeongnam National University of Science and Technology Grant in 2018–2019. BJK was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2017R1A2B2005957). We appreciate the anonymous referee for pointing out that the conventional gravity model plays the role of upper or lower bounds.