A Framework for Spatiotemporal Analysis of Regional Economic Agglomeration Patterns

: Understanding regional economic agglomeration patterns is critical for sustainable economic development, regional resources. Taking Guangdong Province of China as the study area, this paper introduces a comprehensive research framework for analyzing regional economic agglomeration patterns and understanding their spatiotemporal characteristics. First, convergence and autocorrelation methods are applied to understand the economic spatial patterns. Then, the intercity spatial interaction model (ISIM) is proposed to measure the strength of interplay among cities, and social network analysis (SNA) based on the ISIM is utilized, which is designed to reveal the network characteristics of economic agglomerations. Finally, we perform a spatial panel data analysis to comprehensively interpret the inﬂuences of regional economic agglomerations. The results indicate that from 2001 to 2016, the economy in Guangdong showed a double-core/peripheral pattern of convergence, with strengthened intercity interactions. The strength and external spillover effects of Guangzhou and Shenzhen enhanced, while Foshan and Dongguan had relatively strong absorptive abilities. Moreover, expanding regional communication and cooperation is key to enhancing vigorous economic agglomerations and regional network ties in Guangdong by spatial panel data analysis. Our results show that this is a suitable method of reﬂecting regional economic agglomeration process and its spatiotemporal pattern.


Introduction
The spatial agglomeration of economic activities is a typical fact in the process of world economic development. From a global perspective, global production is concentrated in a few large regions. The gross domestic product (GDP) of eight major industrialized countries, including the United States, Japan, Germany, France, the United Kingdom, Italy, Canada and Russia, accounted for more than 57% of the world's GDP in 2017 [1]. At the national level, population and resources are gradually concentrated in cities and urban agglomerations (mega-cities). Projections indicate that by the middle of this century, 66% of the world population would reside in urban areas [2]. In particular, the metropolitan region around Paris has created 29% of GDP in France, using only 2.2% of the total area of France, while the GDP of London and New York metropolitan areas exceeds the GDP of India [3]. In China, the land area of Yangtze River Delta, Pearl River Delta and Beijing-Tianjin-Hebei region is less than 4% of the country, but it gathered 15% of the nationwide population and produced nearly 35% of the country's GDP by 2011 [4]. These empirical facts show that agglomeration is a common phenomenon in the process of economic globalization and regional economic integration.
Driven by the central government's policy or the discovery and utilization of natural resources, such as coal and iron, economic prosperity in the urbanization process is unequally distributed within countries. In these prosperous regions, the sustained expansion of a nonequilibrium economy, social development among regions and the exchange of labor, capital and other resources can be promoted among each other. These resources are prerequisites for regional urban growth and economic agglomeration and form the spatial structural and organic combinations of the economy, society and culture among individual cities [5]. Furthermore, the regional economic core of economic activity communications produces a spillover effect. This can promote the strengthening of the spatial interaction of economic factors in other areas of the region and the upward development of the economy, and cause the region to keep evolving and form urban agglomeration areas [6]. Generally, the regional economy of the urban agglomeration is increasing rapidly and exhibiting more network-like features. However, the over-agglomeration may result in problems, such as deterioration of regional ecological environment, decrease of resources and environmental carrying capacity, and traffic congestion, which can reduce the economic efficiency of the region [7,8]. Quantitatively characterizing economic agglomeration is essential for understanding exchanges and cooperation among areas and the evolution of the regional economy.
To characterize the economic agglomeration in one region, it is important to understand the regional economic agglomeration process and how the agglomeration affects the regional economic activities. As stated by Fujita and Thisse [7], economic agglomeration refers to the effects caused by the agglomeration of social economic activities and related factors in the geographical space. The geographical space can be divided into multiple scales in one region. Strong regional economic differences within the same country imply the existence of agglomerations at the regional scale. Regional agglomeration is also reflected in large varieties of cities, as shown by the stability of the urban hierarchy and the interplay among cities. At a more detailed level, agglomeration arises under the form of various elements, such as traffic and science technology investment in the inner city itself. Economic agglomerations in different scales present diverse features. The spatial distribution characteristics of regional economic agglomeration are prominent at the regional scale, while the main feature of agglomeration is the interplay among cities. In addition, temporality is another important feature because economic agglomeration is a dynamic evolution process. This paper aims to provide a comprehensive method for analyzing regional economic agglomeration patterns and understanding their spatiotemporal characteristics. To address the problem of spatiotemporal interplay among cities, this study proposes a novel gravitation model named the intercity spatial interaction model (ISIM). A new analysis framework involving convergence, autocorrelation analysis, social network analysis (SNA) based on the ISIM, and spatial panel data model is designed to characterize the regional economic spatial pattern and its pattern of spatiotemporal evolution at different scales. This research can provide a reference for government departments and the public to understand economic spatial patterns and the direction of their development, which is expected to give a basis for regional economic integration development decisions.
The remainder of this paper is structured as follows. The Section 2 highlights the results of previous related studies. The Section 3 describes the research area and data. The proposed framework and critic research methods are presented in the Section 4. The Section 5 provides the main findings, including global analysis and networking structure analysis, as well as the output of the estimated panel data regression models. The Section 6 sums up the results and makes a discussion. Section 7 concludes the study and provides the possibilities for future research.

Literature Review
As one of basic elements of cities, economics should be necessarily addressed in order to attain sustainable development [9,10]. Issued by The New Urbanization Plan (NUP, 2014-2020) strategy introduced in March 2014, China emphasizes economy development is an essential prerequisite of regional sustainable development [11]. Many studies have been conducted to test the impact of agglomeration on urban and regional economics and show that the existence of an inverted-U relationship between agglomeration and economic [8,[12][13][14]. In the early stage of economic development, a large number of production factors flow into some cities, and capital, energy, resources, labor, and science and technology continue to gather in some certain cities, named core cities. It has an agglomeration effect, which promotes the leaping development of the economy of the region and can directly promote the rapid growth of the regional and national economy. However, when the core cities have developed to the limits of development that the regional economics can supported, a large number of production factors from non-core development cities will still continue to be absorbed and transmitted to core development cities if there is no balance policy. At this time, it could lead to a decline in the region economy. This is also a situation that the whole society is worried about, that is, the "Matthew effect": core cities are more and more developed, and backward cities are becoming more and more backward, eventually leading to social instability and economic collapse [15]. Additionally, economic agglomeration can cause the redistribution of production factors and resources involving land, which are the targets of urban planning [16]. Therefore, understanding regional economic agglomeration patterns is critical for sustainable economic development, urban planning and proper utilization of regional resources.
As the most striking feature of the geography of economic activities, regional economic agglomeration has always been an important topic in economic geography research. Scholars mainly focus on the nonhomogeneous economic spatial distribution patterns of raw materials, labor, techniques, products, and influencing factors and the processes of agglomeration and spillover effects in the study region [17]. This study began in the 1920s; earlier studies on regional spatial economic agglomerations mainly focused on the fundamental qualitative analysis of unbalanced evolutionary processes and spatial interactions between pole centers and separate cities. Along with theoretical underpinnings (such as the urban-system interaction theory of Zipf [18]; Reilly's gravity model [19], which is applied to retail relations and the central place theory; the spatial interaction theory raised by Ullman and Boyce [20]; and the Krugman's core-peripheral model [21]), studies on the related factors of economic spatial agglomerations have been developed from a single factor function to a multifactor union function. Krugman emphasized that transportation provides a formation condition for economic agglomeration, and different transportation costs lead to different centripetal force and geographical interactions [22]. Feldman pointed out the prominent roles that technological changes and knowledge spillovers play in spatial economic interactions and agglomerations [23]. Quigley found that the acceleration of urban economies is indispensable for the rapid increase of city populations and urbanization [24].
Over the last decade, the application of quantitative models with a spatial dimension has gained prominence. Scholars have mainly focused on the empirical spatial interaction of attribute data in neighboring regions, which means that based on the autocorrelation characteristic of spatial data, the spatial economic structure and related factors are analyzed for discrete data of different geographical units or temporal variations in geographical units. Research perspectives have included the following three aspects: (1) the combination of improved conventional methods and exploratory spatial data analysis (ESDA) methods to measure regional economic interaction mechanisms and variation degrees, where the gravity model, economic membership degree, accessibility and distance decay are improved to measure the spatial economic interaction strength and functional connection among urban areas at different scales; (2) the exploration of determinants by considering economic agglomeration, the spatial vector autoregression model (VAR) and spatial economic panel data models, including the spatial lag model (SLM), spatial error model (SEM) and spatial Durbin model (SDM), which are constructed to investigate the spillover effects among regions under the consideration of a spatial autocorrelation test and time variations in the regression analysis of influencing factors [25,26]; (3) given the spatial ties among economic agglomerations, research has been developed from single-city-based economic ties to multicity-based economic ties. Studies on economic networks mainly include regional economic explicit networks and evolutionary mechanisms in relation to the traffic, trade, information and technology of a region [27,28]. Social network analysis has been applied to discover the inner relationship of economic agglomerations in recent years [29,30].
However, traditional studies oversimplify the complexity of the regional economic agglomeration process by neglecting certain undefined factors of influence in spatiotemporal environments. Most measures of the economic interplay among cities usually consider travel distance and GDP and neglect other factors. Other factors, including the populations and scales of cities, are related to urbanization and urban sprawl, which can give a positive or negative impact on regional economic agglomerations. In addition, the analysis of regional economic agglomeration usually employs the analysis of static network patterns, despite the fact that regional economic agglomeration patterns are the result of dynamic spatiotemporal processes. Therefore, we consider multifactor interactions to describe the spatiotemporal evolution of regional economic agglomerations.

Study Area and Data
Guangdong Province lies along the southeast coast of China, which is mountainous to the north, and in the south, there are coastal plains, hills and tablelands. It is one of the locations in China's Special Zones and one of the core regions in the Guangdong-Hong Kong-Macao Greater Bay Area. As a coastal province bordering Hong Kong and Macao, Guangdong has marked advantages in utilizing overseas capital and technology via foreign economics and technical cooperation. It is one of the fastest-growing economies and most vibrant provinces in China. The total GDP of the whole province reached over 8 trillion yuan in 2016. Moreover, it has already entered the stage of rapid urbanization, and the urbanization rate has reached 67.76%, ranking first in the country. Meanwhile, the differences in economic development between cities have increased. Among them, the urbanization rate of the Pearl River Delta region in south-central Guangdong Province has reached 83.84%, which is equivalent to the levels of medium developed countries, and has entered the mature stage of urbanization. A report from the World Bank released in 2015 showed that the Pearl River Delta region surpassed Tokyo, Japan and possessed the world's largest population and fastest-growing urban agglomeration [31].
This paper focuses on the cities in Guangdong Province (Figure 1), and the researchers have collected data from 21 cities in Guangdong Province from 2003 to 2012. The data sources were mainly economic data and a vectorized province-zoning map. The economic data mainly came from the Statistical Yearbook of Guangdong and the China Statistical Yearbook from 2004 to 2013, where an important feature was the availability of information not only for per capita GDP but also for a broad array of explanatory variables. The urban zoning map of Guangdong Province was obtained from the 1:40-million terrain database via the national basic geography information center. The economic attribute data were projected on these maps and analyzed by data processing software, such as ArcGIS, Geoda, UCINET and Matlab.

Methods
As illustrated in Figure 2, a three-tier research framework was proposed to characterize the regional economic agglomeration comprehensively. Firstly, the data analyzed in this study included a topographic map of the study area, and the statistical yearbook data should be collected. After data collection, the yearbook data, such as GDP, per capita GDP, traffic mileage and fund investment of research and development (R&D) activity, were extracted into the topographic map. Secondly, from the perspective of region, we applied economic convergence and spatial autocorrelation methods to detect global economic agglomeration characteristics. At the city scale, the ISIM was proposed to reveal the economic networking interaction power, and a social network analysis based on ISIM was developed to analyze the economic networking interaction structure. Additionally, from the perspective of inner city, economic influence elements and their spillover effects on neighboring regions through networks were explored.

Methods
As illustrated in Figure 2, a three-tier research framework was proposed to characterize the regional economic agglomeration comprehensively. Firstly, the data analyzed in this study included a topographic map of the study area, and the statistical yearbook data should be collected. After data collection, the yearbook data, such as GDP, per capita GDP, traffic mileage and fund investment of research and development (R&D) activity, were extracted into the topographic map. Secondly, from the perspective of region, we applied economic convergence and spatial autocorrelation methods to detect global economic agglomeration characteristics. At the city scale, the ISIM was proposed to reveal the economic networking interaction power, and a social network analysis based on ISIM was developed to analyze the economic networking interaction structure. Additionally, from the perspective of inner city, economic influence elements and their spillover effects on neighboring regions through networks were explored.

Regional Scale
Per capita GDP Regional agglomeration characteristics Analyze spatiotemporal pattern of regional economic agglomeration Economic convergence analysis

Spatial autocorrelation analysis
Economic spillover effect analysis

Inner City Scale
Economic influence analysis Spatial panel data model

City Scale
Social network analysis

Economic Growth Convergence
Convergence has been the central focus of recent works on regional economic growth. The question of whether poorer regions converge on richer regions has received substantial attention from regional economic growth theories and administrators. The idea of convergence in economics is a popular theory that explains dynamic difference changes in the regional economy. In the 1960s, the concept of steady-state equilibrium was proposed by Solow and Swan, who represented a new

Economic Growth Convergence
Convergence has been the central focus of recent works on regional economic growth. The question of whether poorer regions converge on richer regions has received substantial attention from regional economic growth theories and administrators. The idea of convergence in economics is a popular theory that explains dynamic difference changes in the regional economy. In the 1960s, the concept of steady-state equilibrium was proposed by Solow and Swan, who represented a new classical economic growth theory. They assumed that due to capital-diminishing marginal returns, the growth process eventually reaches a steady state, where the growth rate of the per capita output, capital stock and consumption are equal to the exogenously given rate of technological progress. Hence, in a steady state of economic growth to a certain areal extent, the economic growth speed of poor economics tends to be higher than that of richer regions with similar economic conditions. This phenomenon is called economic growth convergence [32][33][34]. The economic growth convergence model includes sigma convergence, beta convergence and club convergence [35]. Sigma convergence is based on a cross-sectional data analysis, which focuses on the discrete status of per capita income for different economic entities. If the degree of per capita income among various economic entities declines gradually over time, sigma convergence is considered to have occurred. The sigma convergence was described as follows: where y it represents the per capita GDP of the economic unit i at the moment t in the region, and σ t represents the convergence index of the region at the time of t. If σ T+t < σ T appears t years after the year T, it means that this region has regional economic convergence.

Global Spatial Autocorrelation Analysis
Global spatial autocorrelation analysis is used to measure the correlation degree of observations in nearby locations and detect the existence of clustering patterns in space [36]. Moran's I is a primary measure of global spatial autocorrelation, with values ranging from −1 to +1. To identify the high-value clusters of per capita GDP, Moran's I was used to estimate the agglomeration degree. The formula for Moran's I was written as follows [37]: where N is the number of spatial units indexed by i and j; x is the variable of interest; x represents the mean value of x; w ij is a matrix of spatial weights with zeros on the diagonal (w ii = 0); and W is the sum of all w ij .

Intercity Spatial Interaction Model
In this paper, we introduced a modified gravity model to quantify the interaction strength of urban economies due to its similarity with universal gravitation rules and the distance attenuation effects of interactions [38]. In 1931, William J. Reilly first introduced the gravity model for spatial interaction problems in regional economic research. Haynes and Fotheringham proposed the conventional gravity model in 1984 to measure the trade scale between two regions by following the rule of decay with distance. Subsequently, after many amendments and improvements, the model has been widely applied as the most common formulation for measuring bilateral relations.
To generalize, let the scale of each city be represented by M and the distance between cities be represented by d ij . We divided each pair of cities into an origin and a destination, which are notated using the subscripts i and j, respectively. In the original gravity model, the interaction between any pair of cities (F ij ) can be defined as a ratio of the scales of the multiplied cities to the distance between the pair of cities, which was shown as following: where k is a scaling constant, and b represents the friction of the distance parameter controlling the effect of generalized travel time.
To describe intercity spatial interactions precisely, three fundamental modifications need to be made based on the general form of the gravity model. First, Newton's gravitational model presents that the attraction between two heavenly bodies is proportional to the product of their masses. For cities, based on previous studies, the level of economic development and population are important indicators of the mass of a city [39]. In addition, the speed of urban expansion can directly reflect the process of economic development, which is obviously related to changes in the urban land use area within a city [40]. In this paper, to describe changes in urban land use precisely, we used the area of built-up land as one indicator for the size measurement of a city. Thus, GDP, population and a city's size can be applied as important indicators for mass measurements and the economic development level [41,42]. Second, with the development of the economy and infrastructure, the accelerated construction of high-speed railways (HSRs) and highways has brought cities closer together, and the traffic modes of residents have become diverse [43]. The use of travel time can be more accurate for depicting the changes in spatial connections between cities. Therefore, we replaced the distance element (in kilometers) by the travel time (multiplied by 100); by comparing the travel time of highways, railways and HSRs, the minimum travel time was selected in this model. In addition, the friction of distance parameter has proven to be two in the empirical study of Taaffe (1962), which demonstrated that spatial economic interactions are inversely proportional to the square of the distance between two cities. Third, by considering economic imbalance effects, origins and destinations have unequal interaction values proportional to the level of economic development. Hence, the ratio of the origin's GDP to the sum of the GDP for the pair of cities was applied to modify the scaling constant k. Moreover, the relative membership grade (L o ) was applied to divide the level of economic attraction for an origin, and the mean regional interaction strength (M) represents the mean value of the economic interaction.
Therefore, the ISIM can be represented as:

Social Network Analysis Method
From the perspective of relationships and by taking the node tie as an analytical unit, the social network analysis method focuses on the network structure. The network structure depicts the complex attributes of ties between nodes, the diversity of nodes and the inter-relationship between nodes in the network [44]. The SNA is one important method that applies the graph theory and algebra to depict collective interaction patterns and individual behaviors. In this paper, we used the SNA to quantify the influence of the economic association structure of an agglomeration network on regional economic differences. The various cities in a region are considered as points in the network, and the economic spatial interaction between cities is considered to be the tie. These points and lines comprise the regional economic spatial association network [45,46].
The key to characterizing an economic interaction network is gauging the spatial interaction strength and spillover effects of a region effectively and efficiently. The characterizations of economic networking interactions mainly include three aspects: (1) the overall network density, which is applied to characterize the relationships among nodes throughout the whole network, where the density value is between 0 and 1, and the closer the value is to 1, the closer the relationship is between nodes; (2) a measure of the structural centrality, which is used to analyze the central positions of the individuals and groups in the network as a whole; and (3) block modeling, which is analyzed by the convergence of iterated correlations (CONCOR) method, where block modeling was first proposed by Boorman and White [47]. In social network analysis, the subgroups of each factor in the network based on the factor's value are determined to separate the closely related subgroups.

Spatial Panel Data Model
The hypothesis of traditional economic approaches is that each region is independent and regional spatial autocorrelation and spatial heterogeneity are neglected, which can lead to biased results and misleading conclusions [17]. Recently, some models have been developed by considering the impacts of spatial dependence, which can be incorporated in three district ways: the SEM, the SLM and the SDM [17,48,49].
The general forms of SLM and SEM are shown as follows: where y represents an n × 1 vector of observations for the dependent variable y; X denotes an n × k matrix of explanatory variables; β is a k + 1 vector associated with the regression parameter X; and ρ is the regression coefficient of the spatial autoregression, which reflects the spatial dependence that is inherent of the sample observation value. W 1 and W 2 comprise the n × n contiguity matrix, and W 1 y represents the spatially lagged dependent variable. λ is the coefficient of the spatial autoregression, and W 2 ε represents the spatial lag of the error, which reflects the spatial interaction in the error regression structure. µ~N (0, Ω) (the error term µ in the normal distribution) is satisfied; the diagonal elements of the error covariance matrix Ω are: If λ = 0 and ρ = 0, Equation (8) is a normal regression equation, and there is no spatial correlation.
(1) If λ = 0 and α = 0, the mixed autoregressive-regressive model is created and considers spatial dependence across the observations on the dependent variable. Hence, the SLM defined by Anselin [17] accounts for spatial autocorrelation, which can also be called the spatial autoregressive model (SAR); the model is shown as follows: Sustainability 2018, 10, 2800 9 of 22 (2) If ρ = 0 and α = 0, the SEM is created, which can also be called the spatial autocorrelation model (SAC) [50]. The SEM considers spatial dependence across the error terms. The model is defined as follows: In terms of the SLM, a large value of ρ implies considerable spatial dependence and strong influences from neighboring regions. In regard to the SEM, the influences from neighboring regions happen only if the growths of neighboring regions are above or below the normal value of growth. The deviation in the growth of neighboring regions from the expected value refers to the large residual in Equation (12). This effect has a direct correlation with the coefficient λ.
Considering the impacts of both spatially lagged dependent variables and spatially lagged independent variables, the SDM is a better spatial model specification. The advantage of this model is that it produces unbiased and inconsistent estimates. The form of the SDM is presented below: where the interaction terms of spatial weight matrix and exogenous variables are designed to calculate spatial spillover effects among the exogenous variables. All other variables are the same as Equations (8) and (9).
To detect spatial dependence and discriminate the optimal model, five specification tests were introduced in this paper, including Moran's I test, the Lagrange multiplier test for errors (LM-ERR test), the Lagrange multiplier test for lag (LM-LAG test), the Wald test, and the likelihood-ratio (LR) test.

Economic Convergence and Spatial Autocorrelation Structure
Since the reform and open-up policies were instituted in China, the economic development of Guangdong Province has been in a leading position throughout the country. However, at the same time, the economic development of Guangdong shows great regional differences. From a spatial perspective, changes in the spatial distribution of the urban economy in Guangdong Province can be detected by applying the geostatistical approach of area-to-point kriging for the distribution of per capita GDP data collected by 101 county-level cities in Guangdong from 2001 to 2016 ( Figure 3). As we can see from the figure, the overall economic strength of Guangdong grew remarkably. From 2001 to 2011, the Pearl River Delta exhibited the strongest economy, while from 2011 to 2016, Guangzhou and Shenzhen gradually became the economic cores of Guangdong. The economic development of northern mountainous, eastern flanking and western flanking cities was significantly less than that of the Pearl River Delta. Moreover, cities in the Pearl River Delta also showed core and peripheral phenomena. The economic development of Huizhou, Zhaoqing and Jiangmen was also far less than those in other cities of the Pearl River Delta. In other parts of Guangdong, the urban economic development of Zhanjiang, Maoming, Shaoguan, Meizhou and Jieyang was defined by subcenters in the western, northern and eastern regions of Guangdong, while economic development was relatively weak in Meizhou, Heyuan and Yunfu.
The analysis of regional economic spatial convergence can reveal the trend in the overall regional development level and differences in the economic development of Guangdong. As seen from Figure 4, Guangdong showed a general trend of economic convergence. The economy of Guangdong Province was in a more divergent state from 2003 to 2005. Since 2005, the economy has gradually converged and stabilized. The differences among cities decreased constantly, regional economic development gradually moved to an equilibrium state, and the growth of the urban economy in the eastern flanking, western flanking and northern mountainous regions was relatively fast, compared to the cities of the Pearl River Delta. In addition, combined with the values of Moran's I, changes in the spatial distribution of per capita GDP were divided into four periods. From 2001 to 2003, spatial agglomeration in Guangdong strengthened, while regional economic differences decreased gradually, which indicated the increasing speed of non-centralized cities in the study area. During the second period, agglomeration reached a higher level in 2007, and differences among cities slowly became larger; Guangdong was in the core development period from 2004 to 2007. From 2008 to 2014, spatial agglomeration and the differences among cities in Guangdong slowly decreased, which meant that Guangdong was in a period of peripheral and heterogeneous development. During the fourth period from 2015 to 2016, spatial agglomeration began growing, and the differences among cities decreased; Guangdong appeared (again) to exhibit the phenomenon of central-promoting peripheral development.
relatively weak in Meizhou, Heyuan and Yunfu.
The analysis of regional economic spatial convergence can reveal the trend in the overall regional development level and differences in the economic development of Guangdong. As seen from Figure  4, Guangdong showed a general trend of economic convergence. The economy of Guangdong Province was in a more divergent state from 2003 to 2005. Since 2005, the economy has gradually converged and stabilized. The differences among cities decreased constantly, regional economic development gradually moved to an equilibrium state, and the growth of the urban economy in the eastern flanking, western flanking and northern mountainous regions was relatively fast, compared to the cities of the Pearl River Delta. In addition, combined with the values of Moran's I, changes in the spatial distribution of per capita GDP were divided into four periods. From 2001 to 2003, spatial agglomeration in Guangdong strengthened, while regional economic differences decreased gradually, which indicated the increasing speed of non-centralized cities in the study area. During the second period, agglomeration reached a higher level in 2007, and differences among cities slowly became larger; Guangdong was in the core development period from 2004 to 2007. From 2008 to 2014, spatial agglomeration and the differences among cities in Guangdong slowly decreased, which meant that Guangdong was in a period of peripheral and heterogeneous development. During the fourth period from 2015 to 2016, spatial agglomeration began growing, and the differences among cities decreased; Guangdong appeared (again) to exhibit the phenomenon of central-promoting peripheral development.

Analysis of the Networking Structure of Spatial Economic Associations
With the continuous development of the economy, a provincial urban system has gradually formed a networking association structure, as there are obvious differences among cities in Guangdong. The industrial division of labor and cooperation among cities in Guangdong should be further strengthened. Thus, in this paper, we studied the network association relationship in Guangdong Province to deeply explore the status of each prefecture-level city in the agglomeration structure with the aid of the SNA. The SNA is a quantitative and intuitive tool to determine the clear

Analysis of the Networking Structure of Spatial Economic Associations
With the continuous development of the economy, a provincial urban system has gradually formed a networking association structure, as there are obvious differences among cities in Guangdong. The industrial division of labor and cooperation among cities in Guangdong should be further strengthened. Thus, in this paper, we studied the network association relationship in

Analysis of the Networking Structure of Spatial Economic Associations
With the continuous development of the economy, a provincial urban system has gradually formed a networking association structure, as there are obvious differences among cities in Guangdong. The industrial division of labor and cooperation among cities in Guangdong should be further strengthened. Thus, in this paper, we studied the network association relationship in Guangdong Province to deeply explore the status of each prefecture-level city in the agglomeration structure with the aid of the SNA. The SNA is a quantitative and intuitive tool to determine the clear process behind network formations [51].
In this paper, we calculated the economic interaction strengths of 21 Guangdong prefecture-level cities based on the ISIM. The economic interaction strength of a city was set as 0 to exclude self-loops (not counting the diagonal values in the economic interaction matrix). Taking these 21 cities as network nodes, we connected them via a digraph according to their economic interaction strengths, which we defined as a dichotomous division with a threshold value to transform the matrix into an adjacency matrix. The threshold satisfied the basic requirement for connecting all nodes on the map, and the interaction strength should greater than the regional mean interaction strength of a certain year. Hence

Network Cohesion Analysis
Network cohesion represents the connectedness in a group or at the entire network level; more ties between nodes yield a tighter structure, which is (presumably) more cohesive. Structural cohesion can be interpreted by various measures, including density and geodesic distance [52]. Network density patterns are used to estimate cohesive interconnected relationships among all nodes within the network. Geodesic distance refers to the shortest distance between two factors.
In terms of this research, we can see that the density of the whole network increased from 0.2929 in 2001 to 0.3119 in 2016, and the distance-based cohesion increased from 0.449 to 0.497 based on the UCINET output in Table 1, which implied that the ties among cities in the network grew stronger during the last 15 years, the network of Guangdong was at a moderate level, and the cohesiveness needed to be further strengthened. In addition, the average distance for this region decreased from 1.866 to 1.618, demonstrating that, on average, every city should pass one to two nodes to reach the connection with the target city, and the communication among cities was relatively smooth.

Structural Centrality Analysis
The structural centrality of a network is used to identify the contribution of structural positions to the importance, influence and prominence of an actor in the network [52], which includes the measurement of point centrality and the overall centralization of a graph. Structural centrality can be expressed be three basic indicators: degree centrality, closeness centrality and betweenness centrality. Degree centrality refers to the number of ties to other nodes from of a given node, and a higher value of degree centrality indicates a more central position in the network. Closeness centrality represents the graph-theoretic distance of a given nodes to all other nodes, which is an inverse measure of centrality. Betweenness centrality measures the number of geodesic paths that pass through a node. In addition, centralization signifies the overall cohesion or integration of the graph. Table 2 shows the calculation results of the indicators on the centrality. The outdegree centralization of the economic network was 15.640% in 2001 and 17.577% in 2016, and the indegree centralization was 6.291% in 2001 and 6.805% in 2016. The centralization of the study region continuously increased but was still at a relatively low level. From the perspective of the centrality of the cities, we can summarize the findings into three points. (1) Guangzhou and Shenzhen, as the major cities in the Pearl River Delta, had larger values of outdegree centrality than other cities, which indicated that Guangzhou and Shenzhen were always the core of the network and had increasing external spillover effects on the surrounding cities. In addition, these two cities are the only cities that had larger values of outdegree centrality than values of indegree centrality in Guangdong. This meant that the external spillover effects in these two cities were greater than the acquired effects, and Guangzhou and Shenzhen were dominant in the study area; (2) in contrast, during the past 15 years, Heyuan and Yangjiang had relatively low values of outdegree and indegree centrality, illustrating that Heyuan and Yangjiang had weak ties with other cities and less external economic interaction; (3) due to the strong influence of surrounding central cities, Dongguan and Foshan always had relatively high values of indegree centrality, which were larger than the corresponding values of outdegree centrality; (4) from the view of cluster coefficients, the overall cluster coefficient increased from 0.675 to 0.698, and the cumulative effect in Guangdong was enhanced.
Dongguan, Zhongshan, Yunfu and Foshan, which had higher cluster coefficients tended to create tightly knit groups characterized by a relatively high density of ties; and (5) for betweenness centrality, Guangzhou, Zhaoqing, Dongguan and Yunfu had larger values than other cities, demonstrating that they had greater impacts on exchanges and cooperation with other cities in Guangdong in 2001. In 2016, Shenzhen, Guangdong, Huizhou, Foshan and Dongguan had more influences.
By comparing the economic network centrality between two points in time, we can see that (1) the gap between Guangzhou and Shenzhen narrowed. The growth of centrality in Shenzhen has obviously been accelerated during the past 15 years (especially betweenness centrality), which meant that the functions of economic exchange and cooperation in Shenzhen have strengthened; (2) Foshan, Dongguan and Zhongshan were always second-tier members of the economic agglomeration in Guangdong. The interactions of Foshan, Zhongshan and Qingyuan with other cities enhanced, while the statuses of Dongguan and Jiangmen declined; and (3) the number of cities that had impacts on exchanges and cooperation with other cities increased from 12 to 13 over the past 15 years. The statuses of Huizhou and Foshan in terms of exchanges and cooperation have improved.

Subgroup Analysis
To explore the network structure of underlying subgroups and the extensive cooperative relationships of these subgroups, we partitioned the convergence economic network using the CONCOR method in UCINET. The CONCOR is a method of hierarchical clustering for relational data and is applied to uncover the relational position of cities within the network based on structural equivalence [53]. Figures 7 and 8 show the results of the CONCOR. We can find that cities were grouped into four clusters at level 2 and eight clusters at level 3 in 2001 and 2016, respectively, and the intrinsic cooperation and similarity within subgroups and intergroup relationships were detected from the dendrogram.
In 2001, for the economies of central cities and their closely associated cities in Guangdong and the Pearl River Delta region, the first subgroup included Guangzhou, Dongguan and Shenzhen, which formed a level-2 subgroup when combined with Foshan. The third and fourth subgroups were composed of the remaining cities in the Pearl River Delta, of which connections with the center were relatively weaker than those with Foshan. Subgroup five was the northeast economic group, and subgroup six was the western economic group. In addition, subgroups seven and eight comprised the eastern economic group. Combined with the density value, we can see that subgroups one, two and three had closer connections with each other. Subgroup one had full coverage in connection with the other city clusters in Guangdong, while the fifth, sixth and seventh subgroups had very slight economic ties with the other groups (particularly the Meizhou-Shanwei subgroup).
In 2016, the correlation between Guangzhou and Shenzhen strengthened and formed a subgroup, where more frequent exchanges and cooperation occurred between their peripheral cities, including Dongguan, Zhongshan and Foshan; these two subgroups comprised the first division of the economic network at level 2 and had increasing ties with other city clusters. The remaining cities in the Pearl River Delta and Yunfu formed the third subgroup, which indicated that Yunfu grew rapidly, and cities in the Pearl River Delta experienced strengthened correlations and stable correlation structures during the past 15 years. Due to the increasing interactions of the north city cluster with the Pearl River Delta, subgroup three and subgroup four (composed of Shaoguan and Qingyuan, respectively) formed a subgroup at level 2. The fifth to seventh subgroups comprised the eastern cities of Guangdong, and the eighth subgroup comprised the western cities of Guangdong. Compared with the city clusters in 2001, the compositions of the city clusters in the east and north changed slightly; this phenomenon mainly occurred in Heyuan, Shanwei and Meizhou, while the clusters of Jieyang, Shantou and Meizhou and the clusters of Yangjiang, Zhanjiang and Maoming had relatively stable correlation structures. Qingyuan, respectively) formed a subgroup at level 2. The fifth to seventh subgroups comprised the eastern cities of Guangdong, and the eighth subgroup comprised the western cities of Guangdong. Compared with the city clusters in 2001, the compositions of the city clusters in the east and north changed slightly; this phenomenon mainly occurred in Heyuan, Shanwei and Meizhou, while the clusters of Jieyang, Shantou and Meizhou and the clusters of Yangjiang, Zhanjiang and Maoming had relatively stable correlation structures.   cluster with the Pearl River Delta, subgroup three and subgroup four (composed of Shaoguan and Qingyuan, respectively) formed a subgroup at level 2. The fifth to seventh subgroups comprised the eastern cities of Guangdong, and the eighth subgroup comprised the western cities of Guangdong. Compared with the city clusters in 2001, the compositions of the city clusters in the east and north changed slightly; this phenomenon mainly occurred in Heyuan, Shanwei and Meizhou, while the clusters of Jieyang, Shantou and Meizhou and the clusters of Yangjiang, Zhanjiang and Maoming had relatively stable correlation structures.

Analysis of Influential Elements in the Spatial Economy
Apart from the analysis of intercity spatial economic interaction relationships from a regional standpoint, in order to precisely understand the influences of related elements on economic development when considering spatial agglomeration and spillover effects from the city itself, we needed to further explore the economics of Guangdong by using a spatial econometrics estimator for the panel data, as suggested by Elhorst [54].
The development of the regional economy is not only related to its economy structure and industrial composition but also other factors and characteristics of the city, including transportation, the construction of urbanized infrastructure, the advancement of science and technology, geographic positions and the economic spillover effects from surrounding cities. Hence, in this paper, as the main indicators of cities, urbanization rate, fund investment of research and development (R&D) activity and traffic mileage were selected as the independent variables for the spatial panel data model [55][56][57]. As urbanization is an integral part of growth and development transition of economics, the urbanization rate can be used to express population concentration along with economic development [58]. R&D expenditure reflects levels of regional innovation cost and knowledge-driven economy [59]. Traffic infrastructure is the fundamental feature of regional economy communications, one of the representative index is traffic mileage [60]. Details of the variables are outlined in Table 3. In addition, the geographic adjacency relation based on the rook method was adopted to construct a spatial weight matrix. When area i adjoins with area j, the spatial weight matrix (W ij ) is equal to 1; otherwise, W ij = 0. Therefore, the spatial panel data model was presented as follows: (17) where URB it is the urbanization rate; L it is the regional traffic mileage; RD it is fund investment of R&D activity; β 1 , β 2 , β 3 , γ 1 , γ 2 , and γ 3 are the parameters of independent variables URB it , L it , and RD it , respectively; and W represents the NT × NT weighted partition matrix. Each sub-block on the diagonal line of the matrix represents a 21 × 21 binary adjacency square matrix, and all sub-blocks in the nondiagonal area are equal to 0, which represents the locational relationship among the economic subunits in the study area. Equations (15) and (16) represent the general form of the SLM and SEM, and Equation (17) represents the SDM of this research. This paper used observation values from 2001 to 2016 from 21 prefecture-level cities in Guangdong and the spatial econometric models above to generate the weight matrix and calculate the regression equation in Matlab. By calculating the ordinary least squares (OLS) estimation of the total annual per capita output value (y) compared to the urbanization rate (URB it ), regional traffic mileage (L it ) and fund investment of R&D activity (RD it ), the results are shown in Table 4, showing that all estimated coefficients had strong significance. According to the goodness of fit (R 2 ) and adjusted goodness of fit (R a 2 ), the OLS regression model explained nearly 68% of the total variation. For the OLS regression, Moran's I index was used to inspect the spatial autocorrelation in the regression error term. The value of Moran's I was 0.42422661, and there was a significant spatial autocorrelation (p > 0.99) in the residual error. Hence, due to the autocorrelation that existed in the spatial data, the traditional OLS regression model did not fit well and caused the estimation value to deviate or become invalid. Therefore, we applied the SEM and spatial autoregressive model. It can be known from the research of Anselin that the applicability of the SEM and SLM can be determined by using inspection methods via the LM-ERR and LM-LAG [17]. If LM-LAG > LM-ERR, the SLM is chosen; otherwise, the SEM is chosen. If the discriminant results cannot be obtained by using the method of a Lagrangian multiplier, the comparison between the robust LM-LAG and robust LM-ERR can be conducted further. In this paper, the results showed that the LM-LAG value was 97.13876711, which was obviously higher than the LM-ERR value (91.27106331), and Robust LM statistics were in favor of SLM. They both had a relatively strong significance (p < 0.01), which meant that both SLM and SEM were not rejected. According to the results in Table 5, the results obtained from the SDM were relatively steady. To test the hypothesis whether the SDM could be simplified to the SEM, H 0 : γ + ρβ = 0 (γ = [γ 1 , γ 2 , γ 3 ] T ; β = [β 1 , β 2 , β 3 ] T ), we performed the Wald test (644.0684, p < 0.0001) and LR test (115.1846, p < 0.0001), which indicated that the null hypothesis was rejected. Similarly, the null hypothesis that the SDM could be simplified to the SLM (H 0 : γ = 0) was also rejected (Wald test: 119.10956928, p < 0.0001; LR test: 101.4815, p < 0.0001). It implied that both the SEM and the SLM must be rejected in favor of the SDM. In addition, all variables were significant at the 1%, 5% or 10% significance levels, except W × L it . L it , URB it , W × URB it , RD it and W × RD it were positive and highly significant. Considering the erroneous conclusions that analysis based on point estimates of the SDM might cause, we therefore turned to analyze the direct, indirect and total effects. The direct and indirect effects can be used for discovering the existence of spatial spillover effects of exogenous variables and finding the implications of the parameters in the SDM. The results in Table 6 showed that the direct, indirect and total effects of highway mileage and R&D expenditures were positive and significant. It reflected that the rise of highway mileage and R&D expenditures would strongly affect both local economy development and its neighboring economy development. As the connected medium of cities, highway network would promote economic common prosperity between cities through material exchanges. The application and popularization of science and technology of a city could promote regional efficiency economic development, and stimulate technological changes in neighborhoods. In view of the urbanization rate, its direct and total effects were positive and significant while its indirect effect was negative. A possible explanation is that urbanization with one's own city has a positive effect on economy development; while when urbanization rate rises in neighboring cities, the neighboring cities gradually attract the population and a city cannot get more immigrants to support economy development, which then drives the urbanization rate of a city to decrease consequently. Therefore, the rise in the investment of urban construction could increase the differences between cities, which may stimulate the formation of regional multilayer structure over the long term. Hence, the rise in the investment of urban construction, traffic network, science and technology would help significantly increase the regional economy, reasonably promoting the changes of the three factors, which are closely related to the healthy regional development.

Discussion
By utilizing the approaches of economic convergence, SNA based on the proposed ISIM and spatial panel data model, we performed an in-depth exploration of the spatiotemporal characteristics of economic agglomerations in Guangdong Province from 2001 to 2016.

1.
From 2001 to 2016, with the rapid development of the economy and urbanization, the development level of Guangdong Province continually increased towards that of developed countries; the economics showed a double-core/peripheral (i.e., the Guangzhou-Shenzhen peripheral pattern of convergence), and the degree of agglomeration maintained a relatively high level and reached its highest value in 2007. Similarly, Liang [61] verified that a rapid growth is the main feature of Guangdong which has still not reached the optimal sustainable state due to the lack of efficiency.

2.
From the perspective of the whole network, moderate economic ties in this network formed, which implied that this region had a complete network with smooth communications and exchanges. From a regional perspective, the radiation strength and effective central attraction effects of Guangzhou and Shenzhen continuously increased, especially for peripheral cities in the Pearl River Delta. The increase in the centrality of Shenzhen obviously accelerated, with an economic strength that was expected to advance that of Guangzhou. Meanwhile, associations among cities within the Pearl River Delta gradually strengthened. However, Heyuan and Yangjiang had weak ties with other cities and fewer external economic interactions. In addition, among the cities in the Pearl River Delta, Foshan and Dongguan had relatively strong absorptive abilities. The central state of Dongguan has declined, and the importance of Foshan and Zhongshan had improved. Huizhou and Foshan played very important betweenness roles, with exchange and cooperation statuses that had improved.

3.
Guangdong could be divided into eight significant subgroups in 2001 and 2016, with obvious heterogeneity and several changes in composition when compared with the subgroups in 2001, particularly in the western, eastern and northern mountainous areas of Guangdong. From the perspective of subgroup interactions, the economic patterns driven by the subgroups of central cities were formed, and exchanges and cooperation among subgroups needed to be denser.

4.
From the perspective of the city itself from 2001 to 2016, regional economic growth had obvious spatial spillover effects, and low-value clusters had been constantly impacted by high-value clusters. The traffic mileage and fund investment of R&D activity had relatively great impacts on economic agglomeration and networking interactions. The urbanization rate contributed to the local economic development, while it had a negative impact on neighboring cities' economy, due to the siphon effect of the large cities. Urban quality improvements, road construction and technological development could make the exchange of commodities and the work force more frequently and constantly facilitate localized spillovers via local networks. Many empirical studies have verified that economic growth is closely correlated with urbanization, the traffic mileage and fund investment of R&D activity. For example, Xu Zhang [62] studies links between technology level and economy growth. Lee [63] examines the relation between urbanization and economic growth in India.
From the research, we can conclude that the regional node-to-node network structure is like a crystalline form in physics. Megacities tend to be at the core of the crystalline structure, playing a leading role in this region. The big cities are at the key positions and the medium cities are the important nodes to connect the small cities in the crystalline structure [64]. Different levels of node cities are connected through the transportation network. In Guangdong, cities in the Pearl River Delta can be the core or key nodes of the region. However, cities in the western, eastern and northern mountainous areas of Guangdong still cannot form the multilayer structure due to the inapparent regional differences. Hence, it is particularly important to develop a number of regional central cities. According to the analysis results, Maoming, Qingyuan and Jieyang have the advantages in states of economy, centrality and better connections with cities in the Pearl River Delta. These cities can be promoted into regional central cities by increasing innovation investments and population migration, which can make the region form a relatively stable crystalline structure.
Therefore, the government should continuously take active measures to promote the steady economic growth of core cities in Guangdong and constantly expand their radiation ranges. By expanding regional communication and cooperation, the broken barriers of subgroups in Guangdong are the second crucial point. It was found that reducing travel time further by traffic infrastructure construction is the prerequisite for the coordinated development between cities' economy. Hence, it is recommended that policymakers should pay more attention to reducing travel time via urban infrastructure construction in the eastern and western cities of Guangdong and enhancing the absorption abilities of these cities. Furthermore, enhancing technological innovation is the third key point for enhancing vigorous economic agglomerations and providing kinetic energy for healthy economic growth in Guangdong.

Conclusions
Sustainability of economy has been paid significant attention to due to the rapid development of rapid differential development of urban economy, which is related with rational allocation of resources, such as transportation infrastructure, technology, etc. This paper established a three-tier research framework of spatiotemporal analysis of regional economic agglomeration patterns. First, the convergence and autocorrelation methods were applied to understand the economic spatial pattern based on the pattern of spatiotemporal evolution at the regional scale. Then, a novel gravitation model named the ISIM was proposed to measure the strength of the interplay between cities and SNA based on the ISIM, which was designed to reveal the characteristics of the inner rationales of economic agglomerations at the city level. Finally, we performed a spatial regression analysis of the spatial panel time series data to comprehensively interpret the influences of regional economic agglomerations at the inner city scale.
Based on the perspectives of geographical proximity and the network effect, this paper investigated the economic spatial interaction rules and agglomeration pattern in Guangdong Province, China. Our results showed that this is a suitable method of reflecting the regional economic agglomeration process and its spatiotemporal pattern by combining spatial autocorrelation and network analysis. Compared to traditional interplay measures among cities, which typically consider only travel distance and GDP, the ISIM proposed in this study integrated multi-factors including populations and scales of cities, GDP and travel time. The developed research framework can quantity the characteristics of regional economic agglomeration feasibly, because it describes regional economic agglomeration from three scales: global level (regional spatial autocorrelation), city (interplay among cities) and inner city level (spatial autocorrelation in cities). Our new methods improve upon traditional evaluation methods, and not only enrich the approaches used in regional economic agglomeration but also make the analysis results more comprehensive and reasonable. Notably, the method improves both the objectivity of region and urban spatial decision-making and the scientific validity of urban development. Moreover, it may have considerable value in guiding the compact and multi-center development of local urbanization systems.
One advantage of this study is to provide an approach framework for exploring the complicated problem of economy agglomeration patterns from a spatiotemporal perspective. The main findings of this research and their associated applications are relatively thought-provoking. They may contribute to re-guiding regional policy-making, sustainable economic development and urban planning in terms of balance policies for economic agglomeration. However, some limitations still exist in the current research framework. The relationships among economic agglomeration at different scales need be analyzed in depth. A follow-up study to analyze the influence mechanism of economic agglomeration at different scales is required. Additionally, there is a need to further study other influence factors in economic agglomeration, such as culture, access to resources, or avoidance of natural disasters. In future studies, factors other than GDP, transportation, population and scale of cities should be considered, and the results would be more profound.