Spatiotemporal evolution and associated factors of coupling coordination between socio-economic development and ecological health: county-level evidence from the Central Yunnan Urban Agglomeration, China

1 Introduction

Urbanization reshapes regional systems through the expansion of built-up areas, the transformation of urban–rural spatial structures, and the redistribution of population, economic activities, and spatial resources (Chen et al., 2020; Rao and Gao, 2022; Yang et al., 2023). Industrial restructuring, infrastructure development, and population mobility can improve economic productivity and public-service provision, but they may also intensify land conversion, resource consumption, and ecological pressure (Wang et al., 2025; Yang et al., 2025). These tensions are particularly evident in urban agglomerations, where neighboring administrative units are increasingly connected through economic, infrastructural, environmental, and social interactions. Large urban agglomerations can generate substantial development spillovers to surrounding areas through enhanced factor mobility, economic linkages, and regional connectivity, while simultaneously producing uneven development outcomes and environmental pressures (Fu et al., 2022; Li et al., 2021; Zhao et al., 2024a). The Central Yunnan Urban Agglomeration, one of China’s nationally prioritized urban agglomerations, is located in a low-latitude plateau region characterized by mountainous terrain, ecological sensitivity, ethnic diversity, and substantial disparities in development capacity (Zhang et al., 2021). Kunming functions as the dominant political and economic core, whereas many surrounding counties have weaker infrastructure, limited fiscal capacity, and stronger ecological constraints. Assessing whether socio-economic development advances in coordination with ecological health is therefore essential for understanding the sustainability of regional urbanization in this mountainous and multi-ethnic context.

Existing studies have provided important insights into the coordination between human activities and ecological or geographical systems, but several limitations remain (Liu et al., 2018; Zeng et al., 2022). Recent studies have extended coupling-coordination research by integrating multidimensional assessment, spatial diagnostics, and interpretable models to reveal substantial spatial and temporal heterogeneity in socio-economic–ecological relationships across urban agglomerations (Chen et al., 2026; Ma et al., 2025). First, much of the literature has been conducted at the provincial, river-basin, or prefecture level, which can obscure substantial differences among county-level units, particularly in mountainous regions where neighboring counties vary considerably in accessibility, administrative capacity, land-use conditions, and ecological sensitivity (Bi et al., 2023; Chen et al., 2025). Second, many coupling-coordination studies focus primarily on index construction and descriptive spatial patterns, while giving less attention to the distinction between temporal changes within individual regions and persistent structural differences across regions. Conventional global models may overlook local spatial variation, whereas purely spatial models cannot by themselves identify temporal dynamics. Third, some associated-factor analyses introduce explanatory variables that are identical or closely related to the indicators used to construct the coupling coordination degree, creating a risk of circular explanation and mechanically inflated associations. Cultural characteristics are also often discussed through theoretical interpretation or map comparison but are less frequently incorporated into formal statistical and predictive models (Dearing et al., 2014; Li et al., 2025). These limitations indicate the need for a framework that combines intertemporally comparable index construction, circularity control, within–between decomposition, spatial diagnostics, and leakage-controlled predictive analysis.

To address these gaps, this study integrates county-level socio-economic statistics, remote-sensing indicators, cultural-resource data, and geospatial information for 49 county-level units in the Central Yunnan Urban Agglomeration from 2010 to 2023. Socio-economic development and ecological health are evaluated using a pooled entropy-weighted TOPSIS framework, and a baseline coupling coordination degree (CCD-B) is used to describe the overall spatiotemporal evolution of socio-economic–ecological coordination. To reduce mechanical overlap between the dependent variable and subsequent explanatory variables, a reconstructed non-overlapping measure (CCD-R) is recalculated and used as the principal outcome in the associated-factor analyses. The study combines spatial autocorrelation analysis, two-way fixed-effects and Mundlak models, OLS–GWR–MGWR comparisons, and leakage-controlled XGBoost–SHAP analysis to pursue three objectives: (1) to measure the spatiotemporal evolution and spatial differentiation of socio-economic development and ecological health; (2) to assess the evolution and spatial distribution of coupling coordination and evaluate its sensitivity to alternative indicator constructions; and (3) to examine the within-county, between-county, spatially varying, and nonlinear associations between CCD-R and selected socio-economic, infrastructure, accessibility, and cultural characteristics. Cultural factors are treated as contextual explanatory dimensions, and the model results are interpreted as statistical associations, predictive contributions, and observed turning points within the sample range.

2 Research area and data sources

2.1 Research area

The Central Yunnan Urban Agglomeration is situated in the central part of Yunnan Province (Figure 1), China (100°43′–104°49′ E, 24°58′–25°09′ N). Its topography is predominantly mountainous, interspersed with intermontane basins. The study area covers the cities of Kunming, Qujing, Yuxi, Chuxiong, and seven northern counties of Honghe Prefecture, encompassing a total of 49 county-level units. The region spans approximately 111,300 km² (accounting for 29.05% of Yunnan’s total area). In 2023, its population reached 21.99 million and its GDP amounted to 1.95 trillion yuan, representing 47.24% and 61.96% of the provincial totals, respectively. Characterized by a high concentration of population and economic activity alongside a fragile ecological environment, the region serves as a typical case for examining the coupling relationship between socio-economic development and ecological health.

2.2 Construction of the indicator system

2.2.1 Socio-economic development evaluation index system

To evaluate socio-economic development across the Central Yunnan Urban Agglomeration, a baseline indicator system comprising eight indicators was constructed. The selected indicators represent economic scale and productivity, population size, industrial structure, social welfare provision, consumption activity, and the intensity of human activity. The baseline system was used to characterize the overall spatiotemporal evolution of socio-economic development and to calculate the baseline coupling coordination degree (CCD-B). The indicator definitions and pooled entropy weights are presented in Table 1. To minimize mechanical overlap between the constructed outcome and the explanatory variables used in the associated-factor analyses, a reconstructed non-overlapping socio-economic index was additionally developed for the calculation of CCD-R, as described in Section 3.1. The entropy weights recalculated under the alternative indicator specifications are reported in Supplementary Table 1.

Table 1. Baseline socio-economic development indicator system and pooled entropy weights.

Entropy weights were calculated from the complete pooled panel of 686 county-year observations and were held constant across 2010–2023. The weights within the indicator system were normalized to sum to one. Indicators marked with † were included in the baseline socio-economic development index used for the descriptive spatiotemporal analysis but were excluded from the reconstructed non-overlapping index used to calculate CCD-R. Regional GDP was also excluded from CCD-R because of its direct relationship with GDP per capita and population size. The entropy weights recalculated under the alternative indicator specifications are reported in Supplementary Table 1.

2.2.2 Ecological health evaluation index system

The ecological health indicator system captures the multidimensional state of the ecosystem across productivity, service capacity, organization, and resilience. Recent ecological assessments have increasingly incorporated ecosystem-service values, human-threat exposure, and spatially explicit ecosystem-service provision processes to provide a more comprehensive understanding of ecological security and ecosystem functioning (Fang et al., 2024; Yin et al., 2023). Within the constraints of consistent annual county-level data availability, the present indicator system focuses primarily on ecosystem productivity, service capacity, landscape organization, and resilience. Based on prior studies, five indicators were selected (see Table 2).

Table 2. Ecological health evaluation index system.

Entropy weights were calculated from the complete pooled panel of 686 county-year observations and were held constant across 2010–2023. The weights within each indicator system were normalized to sum to one.

2.3 Data sources

The data used in this study comprise socio-economic statistics, remote-sensing data, and geospatial datasets for the period 2010–2023. Socio-economic data were obtained from the Yunnan Statistical Yearbook, the China County Statistical Yearbook, the statistical bulletins of national economic and social development for the 49 counties, and the Fifth, Sixth, and Seventh National Population Censuses. Land-use data, ecosystem service capacity, and landscape diversity indices were obtained from the Resource and Environment Science Data Center, Chinese Academy of Sciences. Net primary productivity (NPP) was derived from the MODIS MYD17A3H product, while the normalized difference vegetation index (NDVI) was obtained from the National Tibetan Plateau Data Center. Elevation data were collected from the Geospatial Data Cloud platform. Night-time light data were obtained from the DMSP-OLS and NPP-VIIRS datasets provided by the National Oceanic and Atmospheric Administration (NOAA). All datasets were processed and resampled to ensure spatial and temporal consistency for subsequent analysis (see Table 3 for details).

Table 3. Data sources and descriptions.

3 Methods

3.1 Indicator normalization and entropy-weighted TOPSIS

The socio-economic development and ecological health levels of the 49 county-level units were evaluated using the entropy-weighted Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS) method (Qiao et al., 2026). The entropy-weighting procedure assigns greater weights to indicators with higher information variability, thereby reducing reliance on subjective judgments, while TOPSIS measures the relative proximity of each county-year observation to the positive ideal solution.

The final matched dataset formed a balanced panel comprising 686 county-year observations, covering 49 county-level units over the period 2010–2023. No missing values remained in the indicators used for index construction; therefore, no interpolation or statistical imputation was required. All indicators in the socio-economic development and ecological health systems were specified as positive indicators, with higher values representing more favorable conditions. Each indicator was normalized using the min–max method (Equation 1):

(1) zij=xij−mini(xij)maxi(xij)−mini(xij)

where x_ij denotes the original value of indicator j for county-year observation i, and z_ij is the corresponding normalized value. To ensure intertemporal comparability, the minimum and maximum values were identified from the complete pooled panel rather than separately for each year.

The entropy weights were also calculated jointly from the full 2010–2023 panel and were held constant throughout the study period. First, the proportional contribution of observation i to indicator j was calculated as (Equation 2):

(2) pij=zij∑i=1nzij

The entropy value of indicator j was then calculated as (Equation 3):

(3) ej=−1lnn∑i=1npijlnpij

where n = 686, and p_ij ln p_ij was defined as zero when p_ij = 0. The information utility and final entropy weight of indicator j were obtained as follows (Equations 4, 5):

(4) gj=1−ej

(5) wj=gj∑j=1mgj

where m denotes the number of indicators within the corresponding subsystem. The weights were normalized to sum to one separately for the socio-economic development and ecological health indicator systems. Using pooled and temporally fixed weights ensures that changes in the composite indices reflect changes in the underlying indicator values rather than annual fluctuations in the weighting structure.

After the entropy weights were obtained, the weighted normalized matrix was calculated as vij=wjzij. The positive and negative ideal solutions were defined as vj+=maxi(vij) and vj−=mini(vij), respectively. The Euclidean distances between each observation and the two ideal solutions were calculated as (Equation 6):

(6) Si+=∑j=1m(vij−vj+)2, Si−=∑j=1m(vij−vj−)2

The relative closeness of observation i to the positive ideal solution was then calculated as (Equation 7):

(7) Ui=Si−Si++Si−

where U_i ranges from 0 to 1, with a higher value indicating a higher level of development or ecological health. The procedure was implemented separately for the socio-economic and ecological indicator systems, producing the socio-economic development index U₁ and ecological health index U₂, respectively.

To address the potential circularity caused by overlap between the socio-economic indicators and the explanatory variables used in the subsequent associated-factor analyses, two principal versions of the socio-economic development index were constructed. The baseline index retained all eight socio-economic indicators and was used to describe the overall spatiotemporal evolution of socio-economic development and baseline coupling coordination. A reconstructed non-overlapping index was then developed by excluding regional GDP, GDP per capita, total population, local general public budget revenue, and nighttime light intensity. Regional GDP was additionally removed because it is directly related to GDP per capita and population size. The reconstructed index therefore retained the tertiary-industry proportion, the number of beds in social welfare institutions, and total retail sales of consumer goods.

For the reconstructed specification, indicator normalization, entropy weighting, TOPSIS scores, and the corresponding coupling coordination degree were recalculated from the beginning rather than obtained through a simple adjustment of the baseline index. The ecological health index remained unchanged. The coupling coordination degree derived from the reconstructed socio-economic index, hereafter referred to as CCD-R, was used as the principal dependent variable in the panel-regression, MGWR, and XGBoost–SHAP analyses. Additional leave-one-out and joint-exclusion index specifications were used as sensitivity tests to assess the dependence of the results on indicator composition.

3.2 Spatial autocorrelation

Spatial autocorrelation analysis (Zhang et al., 2026) was employed to examine the spatial agglomeration patterns and local heterogeneity of both socio-economic development and ecological health levels, incorporating both global and local spatial association measures. A contiguity-based spatial weight matrix in Stata was utilized. Based on the local Moran’s I and its significance, spatial units were categorized into four types: high-high clusters, low-low clusters, high-low outliers, and low-high outliers.

3.3 Coupling coordination degree

The coupling coordination degree (CCD) model (Wang and Wang, 2025) was adopted to quantify the interaction intensity and coordination level between socio-economic development and ecological health. This model simultaneously captures the degree of coupling between the two subsystems and their overall coordinated state, thereby identifying synergistic or conflicting relationships within the socio-ecological dimension.

The coupling degree C is first calculated to measure the strength of interaction between the two subsystems (Equation 8):

(8) C=2×[U1×U2(U1+U2)2]1/2

Where U1 and U2 represent the comprehensive evaluation indices of socio-economic development and ecological health, respectively. The value of C ranges from 0 to 1, with higher values indicating stronger interaction between the two subsystems.

However, the coupling degree alone cannot distinguish between high-level synergy and low-level co-decline. To address this limitation, the comprehensive coordination index T is introduced (Equation 9):

(9) T=αU1+βU2

Following previous coupling-coordination studies that regard socio-economic development and ecological health as equally important dimensions of regional sustainable development, the contribution coefficients were set as α = β = 0.5 (An et al., 2022; Wan et al., 2021; Zou et al., 2022). The coupling coordination degree (D) was then calculated using the following formula (Equation 10):

(10) D=C×T

where D ranges from 0 to 1 and serves as a comprehensive measure of the coordination level between the two subsystems. The classification criteria for D, which have been widely adopted in coupling-coordination studies, are presented in Table 4 (An et al., 2022; Wan et al., 2021). To assess the sensitivity of the results to the relative importance assigned to the two subsystems, both CCD-B and CCD-R were recalculated under three contribution-coefficient settings: (α, β) = (0.4, 0.6), (0.5, 0.5), and (0.6, 0.4). The resulting continuous values, annual trajectories, county rankings, coordination classifications, and absolute differences were compared with the equal-weight baseline.

Table 4. Classification criteria for the coupling coordination degree (D).

The row-standardized Queen-contiguity matrix was used as the primary spatial-weight specification because it directly captures adjacency among county-level units. To assess whether the results were sensitive to the definition of spatial proximity, Global Moran’s I was additionally recalculated using row-standardized inverse-distance and k-nearest-neighbour spatial weight matrices. The same permutation-based significance procedure was applied to all specifications, and the results of the alternative matrices are reported in Supplementary Table 5.

3.4 MGWR

Minority population share was derived from the 2010 and 2020 population censuses. Values for 2011–2019 were linearly interpolated between the two census observations, while the 2020 census share was carried forward to 2021–2023. The resulting annualized series was used consistently in the annual MGWR specifications and the corresponding VIF diagnostics. Because the values for non-census years are approximations, the associated coefficients are interpreted cautiously.

Multiscale Geographically Weighted Regression (MGWR) (Shi et al., 2025) was employed to examine the spatially heterogeneous associations between CCD-R and the selected explanatory variables. The MGWR model utilized an adaptive bi-square kernel function and the corrected Akaike Information Criterion (AICc) for bandwidth selection. The spatially varying regression coefficients were then visualized to illustrate the direction and magnitude of the local associations across the study area, thereby revealing regional characteristics of socio-ecological interactions (Equation 11).

(11) yi=β0(ui,vi)+∑k βk(ui,vi)xik+ϵi

where y_i is the dependent variable (coupling coordination degree) at location i; (u_iv_i) are the geographic coordinates; β₀ is the intercept; β_k(u_iv_i) is the spatially varying coefficient for the k-th explanatory variable x_ik; and ϵ_i is the random error term.

3.5 XGBoost & SHAP

To complement the MGWR analysis and provide a global assessment of variable importance along with nonlinear response interpretation, the eXtreme Gradient Boosting (XGBoost) algorithm combined with SHapley Additive exPlanations (SHAP) was employed (Gao et al., 2025). The XGBoost model constructs an additive ensemble of decision trees, where each new tree is fitted to the residuals of the preceding model. The predicted output for a given observation is the sum of all individual tree predictions (Equation 12):

(12) yi^=∑k=1Kfk(xi), fk∈F

where f_k denotes the k-th tree in the ensemble, x_i is the feature vector of observation i, and F represents the space of all possible regression trees. The model is optimized by minimizing a regularized objective function that balances empirical loss and model complexity.

To avoid information leakage caused by a random row-wise split, model training and evaluation were conducted using three complementary validation strategies. First, nested spatial GroupKFold validation used five outer folds, with all observations from the same county assigned exclusively to either the training or test fold. Hyperparameters were selected separately within each outer training sample using four-fold county-grouped cross-validation. Second, expanding-window temporal validation was conducted for 2018–2023, with each test year predicted using only observations from preceding years. Hyperparameter tuning within each temporal training window also employed county-grouped cross-validation. Third, a strict spatiotemporal holdout excluded 10 counties from model training and evaluated their observations from 2020–2023 using models trained on the remaining counties through 2019. Model performance was evaluated using R², root mean squared error (RMSE), and mean absolute error (MAE).

Hyperparameters were optimized using a 20-iteration randomized search that minimized cross-validated RMSE. The search space included 200–800 trees for n-estimators, 2–6 for max-depth, 0.01–0.15 for learning-rate on a log-uniform scale, 1–10 for min-child weight, 0.65–1.00 for both subsample and colsample-bytree, 0–0.25 for gamma, 10⁻⁵–1.0 for reg-alpha on a log-uniform scale, and 0.1–20.0 for reg-lambda on a log-uniform scale. All models used the squared-error regression objective, RMSE as the evaluation metric, the histogram tree-building algorithm, and a fixed random seed of 2026. Hyperparameter selection was conducted exclusively within the corresponding training samples, without using information from the outer test observations.

While XGBoost provides overall feature importance, it does not directly reveal the direction or pattern of individual feature contributions. To address this, SHAP values, derived from cooperative game theory, were calculated to decompose the prediction of each observation into additive contributions of individual features. For a given prediction, the SHAP value of the j-th feature is defined as (Equation 13):

(13) ϕj=∑S⊆N∖{j}|S|!(|N|−|S|−1)!|N|![f(S∪{j})−f(S)]

where N is the set of all features, S is any subset of N excluding feature j, and f(S) is the model prediction using the features in subset S. To ensure that the interpretation was based on genuinely unseen observations, SHAP values were calculated only for the test observations in each outer spatial-validation fold and were subsequently pooled as out-of-fold SHAP values. Nonlinear response patterns were estimated using LOWESS smoothing with a smoothing fraction of 0.45. Observed turning points were defined as zero crossings of the smoothed out-of-fold SHAP curves within the 5th–95th percentile range of each predictor. Their stability was evaluated using 300 county-grouped bootstrap resamples, which preserved the repeated-observation structure within counties. Only turning points successfully identified in at least 60% of the bootstrap samples were reported, together with their 95% bootstrap intervals. These values represent sample-specific changes in predictive contribution rather than causal or universal policy thresholds. SHAP values satisfy the properties of local accuracy, missingness, and consistency, providing a theoretically grounded and interpretable decomposition of model predictions.

3.6 Associated-factor framework and circularity control

In addition to the baseline coupling coordination degree (CCD-B), six alternative outcome variables were constructed to evaluate sensitivity to indicator overlap. Four leave-one-out versions were calculated by separately excluding GDP per capita, total population, local general public budget revenue, and nighttime light intensity from the socio-economic development index. A joint-exclusion version, CCD-A, simultaneously excluded these four overlapping indicators. A stricter reconstructed version, CCD-R, further excluded regional GDP because of its direct relationship with GDP per capita and population size. For each specification, indicator normalization, entropy weighting, TOPSIS evaluation, and coupling coordination calculations were repeated from the beginning. The ecological-health index remained unchanged across these specifications. CCD-R was used as the main outcome for the associated-factor models, while CCD-B, CCD-A, and the leave-one-out outcomes were used for sensitivity analysis.

4 Results

4.1 Socio-economic development

4.1.1 Temporal evolution of socio-economic development

As shown in Figure 2, the socio-economic development index increased substantially across all five prefecture-level regions of the Central Yunnan Urban Agglomeration between 2010 and 2023. The regional mean rose from 0.0734 in 2010 to 0.2360 in 2023, representing an overall increase of approximately 221.3%. Kunming maintained a clear leading position throughout the study period, with its mean index increasing from 0.1173 to 0.3327. Qujing and Honghe gradually formed a secondary development tier, reaching 0.2391 and 0.2212, respectively, by 2023. Yuxi and Chuxiong also exhibited continuous improvement, although their 2023 values remained comparatively lower at 0.1858 and 0.1533.

The trajectories reveal both regional catch-up and persistent spatial polarization. Kunming recorded the largest absolute increase, rising by 0.2154, whereas Honghe and Chuxiong experienced the fastest proportional growth, with their indices increasing by more than 280% from relatively low initial levels. Nevertheless, the gap between the highest- and lowest-performing regions widened from 0.0771 in 2010 to 0.1793 in 2023, indicating that the faster relative growth of peripheral regions was insufficient to close the absolute development disparity with Kunming. After sustained growth through 2021, the regional mean remained nearly unchanged between 2021 and 2022, decreasing slightly from 0.2267 to 0.2263, before recovering to 0.2360 in 2023. This temporary slowdown was accompanied by small declines or stagnation in Kunming, Honghe, and Yuxi, while Qujing and Chuxiong continued to increase. Overall, the temporal evolution of socio-economic development was characterized by broad-based growth within a persistent core–periphery structure.

4.1.2 Spatial evolution of socio-economic development

Figure 3 illustrates the spatial evolution of socio-economic development across the Central Yunnan Urban Agglomeration during 2010–2023. In 2010, most county-level units were concentrated in the lowest development category of 0–0.16, while relatively high values were limited to a small number of central urban districts. This distribution revealed a pronounced core–periphery structure, with the main urban area of Kunming forming the dominant socio-economic core and most surrounding mountainous counties remaining at comparatively low development levels. Between 2012 and 2016, the central high-value area became increasingly consolidated, while several prefecture-level urban centers and their surrounding counties gradually shifted from the lowest category to the 0.17–0.32 and 0.33–0.48 categories.

The spatial expansion of medium- and high-value areas became more evident after 2016. By 2018 and 2020, counties with intermediate development levels were distributed more extensively across the northern, eastern, southern, and western parts of the agglomeration, indicating that socio-economic improvement was no longer confined to the Kunming core. Secondary development nodes associated with the principal urban centers of Qujing, Yuxi, Chuxiong, and northern Honghe became increasingly visible. In 2023, the number of counties in the lowest category had contracted substantially, while the 0.33–0.48 and 0.49–0.64 categories expanded across several subregions. The highest development category remained concentrated in the central urban area, however, demonstrating the persistent primacy of Kunming.

Overall, the spatial evolution of socio-economic development was characterized by simultaneous upgrading, outward diffusion, and continued polarization. Although peripheral counties experienced widespread improvement and several secondary development centers emerged, the highest values remained strongly concentrated in the central core. The regional pattern therefore evolved from a predominantly low-level and highly concentrated structure toward a more differentiated hierarchical system comprising a dominant core, emerging secondary centers, and improving peripheral areas.

4.1.3 Spatial autocorrelation of socio-economic development

The Global Moran’s I statistics reported in Table 5 indicate that socio-economic development exhibited significant positive spatial autocorrelation throughout the study period. Moran’s, I remained positive in every year and ranged from 0.3738 to 0.4512, with all permutation-based pseudo p-values below 0.001. The 14-year mean Moran’s I was 0.4186, confirming that counties with similar socio-economic development levels were consistently clustered in geographical space rather than randomly distributed.

Table 5. Global Moran’s I of socio-economic development under the queen-contiguity matrix, 2010–2023.

The pseudo p-values were obtained through permutation tests based on a row-standardized Queen-contiguity spatial weight matrix.

Temporally, Moran’s I increased from 0.3999 in 2010 to 0.4492 in 2013, indicating an initial strengthening of spatial concentration. After moderate fluctuations between 2014 and 2016, the index reached its maximum value of 0.4512 in 2017. This pattern suggests that the spatial differentiation between the highly developed central core and the less-developed peripheral counties became increasingly consolidated during the first half of the study period. After 2017, Moran’s I generally declined, falling from 0.4512 to 0.3738 by 2023. Nevertheless, the value remained positive and highly significant, indicating that the decline represented a weakening, rather than disappearance, of spatial clustering. The gradual reduction in spatial autocorrelation is consistent with the expansion of intermediate development levels and the emergence of secondary development nodes shown in Figure 3.

The local spatial autocorrelation patterns in Figure 4 provide further evidence of this spatial structure. High–high clusters remained concentrated in the central part of the urban agglomeration in 2010, 2017, and 2023, confirming the persistent role of the Kunming urban core as the principal concentration of high socio-economic development. In contrast, low–low clusters were mainly distributed in the northwestern and southwestern peripheral areas and remained relatively stable across the three representative years. These clusters indicate that low development levels in peripheral counties were not isolated phenomena but formed geographically contiguous areas of relative disadvantage.

A small number of high–low and low–high spatial outliers appeared around the boundaries between the central core and the peripheral areas. Their changing locations between 2010 and 2023 indicate increasing local differentiation within the broader core–periphery structure. High–low outliers may represent emerging local growth centers surrounded by less-developed neighbors, whereas low–high outliers identify counties that had not yet fully benefited from the development advantages of adjacent high-value areas. By 2023, the continued presence of the central high–high cluster, together with the increased occurrence of transitional outliers, suggests that spatial diffusion had strengthened but remained uneven. The spatial-autocorrelation results were robust to alternative definitions of spatial proximity. Under the inverse-distance matrix, Moran’s I remained positive and significant in every year, ranging from 0.063 to 0.082, while the k-nearest-neighbor specification produced values ranging from 0.398 to 0.478. Both alternative specifications reproduced the main temporal pattern identified under the Queen-contiguity matrix: spatial concentration strengthened during the earlier period, peaked around 2017, and subsequently weakened while remaining statistically significant through 2023 (Supplementary Table 5). Overall, the Global and Local Moran’s I results jointly demonstrate a persistent yet gradually weakening pattern of positive spatial dependence in regional socio-economic development.

4.2 Ecological health

4.2.1 Temporal evolution of ecological health

As shown in Figure 5, ecological health in the Central Yunnan Urban Agglomeration remained broadly stable but exhibited pronounced interannual fluctuations between 2010 and 2023. The regional mean increased from 0.3393 in 2010 to a period-high value of 0.3624 in 2017, after which it fluctuated moderately and reached 0.3586 in 2022. A noticeable decline occurred in 2023, when the regional mean fell to 0.3438. The net increase over the entire study period was therefore only approximately 1.34%, indicating that ecological health did not follow a continuous upward trajectory. Instead, its temporal evolution was characterized by moderate improvement during the middle of the study period, followed by renewed fluctuations and a recent decline.

Considerable differences were observed among the five prefecture-level regions. Kunming generally maintained the highest ecological health level from 2010 to 2022, rising from 0.3731 in 2010 to a peak of 0.4001 in 2017 before declining to 0.3707 in 2023. Yuxi displayed the most sustained improvement, increasing from 0.3615 to 0.3838 and becoming the highest-ranked region in 2023. Honghe fluctuated around the regional mean and reached its highest value of 0.3678 in 2022, followed by a marked decline to 0.3340 in 2023. Chuxiong exhibited moderate fluctuations, increasing from 0.3238 in 2010 to a peak of 0.3494 in 2017 and ending at 0.3315 in 2023. Qujing remained the lowest-performing region throughout most of the study period, with its index changing only slightly from 0.2853 in 2010 to 0.2836 in 2023. The regional decline observed in 2023 was shared by Kunming, Qujing, Chuxiong, and Honghe, whereas Yuxi continued to improve. These results indicate substantial regional differentiation and temporal instability, but they do not by themselves identify the external factors responsible for the short-term changes.

4.2.2 Spatial evolution of ecological health

Figure 6 presents the spatial distribution of ecological health across the Central Yunnan Urban Agglomeration from 2010 to 2023. Throughout the study period, most county-level units remained within the intermediate range of 0.29–0.38, indicating that the overall spatial structure of ecological health was relatively stable. Nevertheless, considerable local differentiation was evident. Low-value areas, with ecological health indices below 0.28, were persistently concentrated in the eastern part of the urban agglomeration, forming a relatively continuous ecological low-value zone. By contrast, counties with values above 0.39 were mainly distributed in the central part of the study area and in several scattered peripheral locations.

From 2010 to 2020, a relatively stable high-value nucleus was maintained around the central part of the agglomeration. Several counties in this area reached the 0.49–0.58 or 0.59–0.68 categories, while neighboring counties generally fell within the 0.39–0.48 category. However, the spatial extent of high-value areas did not expand continuously. Instead, their distribution fluctuated across the selected years, reflecting the temporal instability observed in the regional ecological health series. In particular, the highest-value category became less prominent by 2023, while the number of counties in the 0.39–0.48 category increased around the central area.

The 2023 pattern displayed a clearer spatial contrast between the central and eastern parts of the agglomeration. Although several central counties continued to exhibit relatively favorable ecological conditions, the eastern low-value zone remained evident and became more spatially continuous. Meanwhile, most western, northern, and southern counties remained in the intermediate category. Overall, the spatial evolution of ecological health was characterized by broad stability, localized fluctuations, and persistent regional differentiation rather than continuous and uniform improvement. The coexistence of a central relatively high-value area and an eastern low-value belt indicates that ecological disparities remained structurally embedded within the regional landscape.

4.2.3 Spatial autocorrelation of ecological health

The Global Moran’s I results presented in Table 6 demonstrate that ecological health exhibited significant positive spatial autocorrelation throughout the study period. Moran’s I remained positive in all years, ranging from 0.3234 in 2010 to 0.5496 in 2023, and all permutation-based pseudo p-values were below 0.001. The mean Moran’s I over the 14-year period was 0.3869, indicating that counties with similar ecological health levels were consistently more likely to be located near one another than would be expected under a spatially random distribution.

Table 6. Global Moran’s I of ecological health under the queen-contiguity matrix, 2010–2023.

Global Moran’s I was calculated using a row-standardized Queen-contiguity spatial weight matrix. Statistical significance was evaluated using permutation tests, and all reported pseudo p-values were below 0.001.

The temporal evolution of Moran’s I can be divided into two broad stages. From 2010 to 2021, the index fluctuated within a relatively narrow range of approximately 0.32–0.39. After increasing from 0.3234 in 2010 to 0.3737 in 2011, spatial clustering remained broadly stable, with only moderate interannual variation. A marked increase occurred after 2021: Moran’s I rose to 0.4811 in 2022 and further to 0.5496 in 2023. This late-stage increase indicates a substantial strengthening of spatial dependence and polarization in ecological health. Importantly, the higher Moran’s I does not necessarily imply an overall improvement in ecological conditions; rather, it indicates that high- and low-value counties became more spatially concentrated.

The Local Moran’s I maps in Figure 7 further reveal the geographical structure underlying the global pattern. Low–low clusters were persistently concentrated in the eastern part of the urban agglomeration in 2010, 2017, and 2023. Their spatial continuity and temporal stability indicate that relatively low ecological health in the eastern counties constituted a regional cluster rather than a set of isolated local conditions. High–high clusters were mainly located in the central part of the study area. The extent of the central high–high cluster increased between 2010 and 2023, indicating stronger spatial concentration among counties with relatively favorable ecological health.

Low–high outliers appeared near the margins of the central high-value cluster in 2010 and 2017, representing counties with relatively low ecological health surrounded by higher-value neighbors. By 2023, these transitional outliers had largely disappeared, while the high–high cluster became more spatially continuous. High–low outliers were rare or absent across the representative years. Overall, the combination of an expanding central high–high cluster and a persistent eastern low–low cluster explains the pronounced increase in Global Moran’s I after 2021. The ecological-health results were similarly robust to alternative spatial-weight specifications. Moran’s I remained positive and significant under both the inverse-distance and k-nearest-neighbor matrices throughout 2010–2023. The inverse-distance values increased from 0.054 in 2010 to 0.103 in 2023, while the k-nearest-neighbor values increased from 0.281 to 0.542. Both specifications reproduced the relatively stable spatial clustering observed through 2021 and the pronounced strengthening in 2022 and 2023, supporting the conclusion that late-stage ecological polarization was not an artefact of the Queen-contiguity matrix (Supplementary Table 5). These results indicate that the ecological health pattern became more spatially polarized toward the end of the study period, even though the regional average remained broadly stable and declined moderately in 2023.

4.3 Coupling coordination patterns and reconstruction sensitivity

4.3.1 Baseline spatiotemporal evolution of coupling coordination

Figure 8 shows a marked improvement in the coupling coordination between socio-economic development and ecological health from 2010 to 2023. Table 7 further quantifies this transition by reporting the number and proportion of counties in each coordination class for 2010, 2014, 2018, and 2023. The mean CCD-B increased from 0.378 in 2010 to 0.516 in 2023. In 2010, 37 of the 49 counties, accounting for 75.51% of the study area, were classified as either severely or mildly imbalanced. The number of imbalanced counties subsequently declined continuously: severe imbalance disappeared after 2010, while mild imbalance fell from 36 counties in 2010 to four in 2018, one in 2020, and none after 2020. Meanwhile, the number of moderately or highly coordinated counties increased from five in 2010 to 19 in 2023. Nevertheless, barely coordinated counties remained the dominant category in the later period, accounting for 63.27% in 2018 and 61.22% in 2023, indicating that the regional system had broadly moved beyond imbalance but had not yet achieved uniformly high-level coordination. The complete annual distribution of counties across the five coordination classes from 2010 to 2023 is reported in Supplementary Table 3.

Table 7. Distribution of counties across CCD-B coordination classes in selected years.

Values are the number of counties, with percentages in parentheses. Each annual cross-section contains 49 counties.

Spatially, the distribution exhibited a persistent core–periphery structure. In 2010, the four highest-ranked counties were Panlong, Wuhua, Guandu, and Xishan districts in Kunming, followed by Hongta District in Yuxi. By 2023, Guandu, Wuhua, Xishan, Panlong, and Chenggong occupied the top five positions, while Hongta District, Anning City, and Qilin District also entered the highly coordinated group. These changes suggest that the high-coordination core remained centered on Kunming’s main urban districts but gradually extended towards adjacent urban centers. In contrast, lower-ranked counties remained concentrated mainly in the peripheral parts of Chuxiong Prefecture, together with several outer counties of Qujing and Kunming, demonstrating that the improvement in coordination did not completely eliminate the regional development gap. Detailed county-level CCD-B and CCD-R values, coordination classes, and rankings for 2010 and 2023 are provided in Supplementary Table 2.

The reconstructed non-overlapping measure, CCD-R, provides a more conservative assessment but preserves the broad spatial hierarchy. Its regional mean increased from 0.374 in 2010 to 0.496 in 2023, and the four highest-ranked counties in 2023 remained Guandu, Wuhua, Xishan, and Panlong districts. Supplementary Figure 1 and Supplementary Table 2 further show that the central high-value structure and the overall upward trajectory were broadly preserved under the reconstructed non-overlapping measure. However, differences in individual county ranks and coordination classes also indicate that absolute classifications are sensitive to the construction of the socio-economic subsystem.

4.3.2 Sensitivity of coupling coordination to indicator reconstruction

The sensitivity of the coupling coordination measurement was assessed by sequentially and jointly excluding socio-economic indicators that overlapped with the explanatory variables used in the subsequent associated-factor analyses. As shown in Table 8, the leave-one-indicator-out specifications remained highly consistent with the baseline CCD-B. Pearson correlations ranged from 0.9948 to 0.9993, Spearman correlations ranged from 0.9855 to 0.9984, and classification agreement remained between 92.27% and 97.23%. The mean absolute differences were all below 0.008, indicating that the exclusion of any single overlapping indicator had only a limited effect on county-level values, rankings, and coordination categories.

Table 8. Sensitivity of the coupling coordination degree to the exclusion of overlapping socio-economic indicators.

CCD-B is based on the original eight-indicator socio-economic development index. CCD-A excludes GDP per capita, total population, local general public budget revenue, and nighttime light intensity. CCD-R additionally excludes regional GDP. All socio-economic indices, entropy weights, and coupling coordination values were recalculated using the complete pooled panel. Classification agreement denotes the proportion of county-year observations assigned to the same coordination category as under CCD-B.

Larger differences emerged under the joint-exclusion specifications, although the broad spatiotemporal pattern remained stable. CCD-A retained Pearson and Spearman correlations of 0.9788 and 0.9655 with CCD-B, while the corresponding correlations for CCD-R were 0.9571 and 0.9254. Their classification agreement rates were 78.57% and 78.72%, respectively, showing that the stricter reconstruction affected absolute categories more strongly than relative ordering. Nevertheless, the annual means of CCD-B, CCD-A, and CCD-R all increased substantially between 2010 and 2023, and the central high-value structure remained visible under CCD-R (Supplementary Figure 1, Supplementary Table 2). These results indicate that the overall improvement in socio-economic–ecological coordination was not solely generated by overlapping indicators. CCD-R was therefore adopted as the principal dependent variable in the subsequent associated-factor analyses, while CCD-B was retained as the baseline measure for describing the long-term spatiotemporal evolution. The spatial comparison further shows that the central high-value structure and peripheral lower-value pattern were broadly preserved under CCD-R (Supplementary Figure 1).

Additional subsystem-weight sensitivity tests showed that the temporal trajectories and county-level rankings of both CCD-B and CCD-R remained nearly unchanged under (α, β) = (0.4, 0.6) and (0.6, 0.4), with all Pearson correlations exceeding 0.998 and annual-trajectory correlations exceeding 0.999. Classification agreement was lower, ranging from 78.13% to 82.94%, indicating that discrete coordination classes were more sensitive than continuous values and rankings to the relative subsystem weights (Supplementary Table 4).

4.4 Associated factors

4.4.1 Panel regression

To examine the temporal and cross-sectional associations between the selected factors and coupling coordination, three complementary panel specifications were estimated using the reconstructed non-overlapping coupling coordination degree (CCD-R) as the dependent variable. The contemporaneous and one-year-lagged models employed county and year fixed effects, with standard errors clustered at the county level. Because distance to Kunming was treated as a time-invariant county characteristic, its coefficient could not be identified in the county fixed-effects models. A correlated random-effects model based on the Mundlak specification was therefore further estimated to distinguish within-county temporal associations from between-county differences and to incorporate the time-invariant distance variable.

The contemporaneous two-way fixed-effects results are reported in Model 1 of Table 9. Road-network density was positively but not significantly associated with CCD-R (β = 0.0054, p = 0.416). In contrast, per capita savings showed a significant negative within-county association with CCD-R (β = −0.0241, p = 0.042). Because the variables were logarithmically transformed and standardized before estimation, this result indicates that a one-standard-deviation increase in per capita savings relative to a county’s own temporal trajectory was associated with a 0.0241 decrease in CCD-R, after controlling for county-specific characteristics and common year effects. The clustered robust joint test of the two explanatory variables was not statistically significant at the conventional level (p = 0.126).

Table 9. Panel estimates and within–between decomposition of factors associated with CCD-R.

Clustered standard errors at the county level are reported in parentheses. All continuous explanatory variables were logarithmically transformed and standardized. Within-county terms represent deviations from each county’s long-term mean, whereas between-county terms represent county-specific temporal means. Distance to Kunming was calculated as the median county-level distance over 2010–2023 and treated as time invariant. *p<0.10; p< 0.05; *p< 0.01.

The one-year-lagged specification produced a consistent pattern. Lagged road-network density remained statistically insignificant (β = 0.0056, p = 0.446), whereas lagged per capita savings retained a significant negative association with current CCD-R (β = −0.0207, p = 0.019). The joint clustered test was marginally significant at the 10% level (p = 0.063). These findings indicate that annual changes in road infrastructure were not systematically associated with contemporaneous or subsequent changes in coordination, while increases in savings were associated with lower CCD-R in both the contemporaneous and lagged specifications. However, these estimates represent within-county statistical associations and should not be interpreted causally.

The Mundlak model further revealed a clear distinction between within-county and between-county relationships. The within-county coefficients were nearly identical to those obtained from the fixed-effects model: road-network density remained insignificant (β = 0.0054, p = 0.400), while per capita savings showed a significant negative association (β = −0.0241, p = 0.035). In contrast, the county-level long-term means of both variables were positively associated with CCD-R. Counties with higher average road-network density exhibited significantly higher coordination levels (β = 0.1027, p = 0.002), and counties with higher average per capita savings also showed higher CCD-R values (β = 0.0536, p = 0.023). This contrast suggests that short-term within-county changes and persistent between-county development differences should not be interpreted as equivalent relationships.

Distance to Kunming was positively associated with CCD-R in the Mundlak model (β = 0.0442, p < 0.001). This result indicates that, after controlling for long-term road density, savings, and common temporal changes, counties located farther from the provincial capital tended to exhibit higher reconstructed coordination levels. This association should not be interpreted as evidence that remoteness improves socio-ecological coordination. Rather, it indicates that the spatial distribution of CCD-R may also reflect persistent characteristics of peripheral counties that are not fully captured by the included variables. The Mundlak model explained approximately 79.1% of the observed variation and was jointly significant at p < 0.001.

4.4.2 MGWR diagnostics and spatial heterogeneity

Before estimating the MGWR models, annual variance inflation factors (VIFs) were recalculated using the same eight explanatory variables included in each annual MGWR specification. As shown in Table 10, GDP per capita, local budget revenue, distance to Kunming, ICH density, and minority population share remained below a VIF of 5 throughout 2010–2023. Nighttime light intensity showed moderate multicollinearity in four years, with a maximum VIF of 7.232, while per capita savings exceeded 5 in seven years and reached a maximum of 6.865. Road-network density exhibited the highest collinearity, exceeding 5 in 12 of the 14 years and reaching 10.083 in 2023. Thus, severe multicollinearity was limited to road-network density in the 2023 cross-section, while the remaining variables remained below the conventional threshold of 10. Minority population share showed consistently low VIF values, ranging from 1.228 to 1.714. Accordingly, all variables were retained to preserve consistency across the annual MGWR specifications, but the road-network-density coefficient in 2023, together with coefficients estimated in years showing moderate collinearity, was interpreted cautiously.

Table 10. VIF diagnostics for the annual core MGWR models and the 2020 cultural-extension model.

VIFs were calculated separately for each annual cross-section using the same eight explanatory variables included in the corresponding annual MGWR specifications. Minority population share was derived from the 2010 and 2020 population censuses, linearly interpolated for 2011–2019, and held at the 2020 census value for 2021–2023. Per capita savings was harmonized to CNY 10,000 per person before logarithmic transformation, while distance to Kunming was treated as a time-invariant county characteristic. VIF values below 5 were regarded as acceptable, values between 5 and 10 indicated moderate multicollinearity requiring cautious interpretation, and values above 10 indicated severe multicollinearity.

To evaluate the suitability of spatially varying regression, OLS, GWR, and MGWR models were compared for 2010, 2014, 2018, and 2023 using the same dependent variable and explanatory-variable specifications. As shown in Table 11, all three models achieved relatively high explanatory power, with adjusted R² values ranging from 0.748 to 0.868. GWR produced the highest adjusted R² in all four years, reaching 0.799, 0.868, 0.842, and 0.765 in 2010, 2014, 2018, and 2023, respectively. MGWR also provided comparatively strong explanatory performance, with adjusted R² values of 0.788, 0.837, 0.816, and 0.761.

Table 11. Comparison of OLS, GWR, and MGWR model diagnostics.

The models were estimated using 49 county-level observations in each year. GWR and MGWR employed an adaptive bi-square kernel. Statistical significance of residual Moran’s I was evaluated using 999 permutations and a four-nearest-neighbor spatial weight matrix.

However, the information criteria did not indicate a consistent advantage of MGWR over the simpler alternatives. Based on the package-reported AICc values, OLS provided the most parsimonious specification in 2010, 2018, and 2023, whereas GWR performed best in 2014. MGWR did not produce the lowest AICc in any of the four years. Residual Moran’s I values were statistically insignificant for OLS, GWR, and MGWR across all selected years, with pseudo p-values exceeding 0.20, indicating that the residual spatial dependence had been adequately controlled even by the global models. Accordingly, MGWR was retained primarily as a supplementary tool for identifying variable-specific spatial scales and local coefficient variation, rather than being interpreted as an unequivocally superior model.

The variable-specific bandwidths further indicate that most associations operated at relatively broad spatial scales. The intercept, GDP per capita, nighttime light intensity, road-network density, per capita savings, and distance to Kunming generally selected bandwidths between 45 and 47 nearest neighbors, accounting for more than 90% of the 49 county-level observations. Their estimated relationships therefore approached global or regional-scale associations, with comparatively limited local variation.

More localized patterns were identified for local budget revenue and minority population share. The bandwidth of local budget revenue was 31 in 2010 and 32 in 2014, indicating stronger spatial heterogeneity during the earlier study period, before increasing to 46 in 2018 and 44 in 2023. Minority population share selected bandwidths of 34, 47, 41, and 45 in the four respective years, suggesting that its spatial association was more localized in 2010 and 2018. Intangible cultural heritage density selected a moderately local bandwidth of 41 in 2014 but approached a near-global scale in the other years (Table 12). Overall, the bandwidth results indicate that strong multiscale variation was concentrated in a limited number of fiscal and cultural variables rather than occurring uniformly across all explanatory factors.

Table 12. Variable-specific MGWR bandwidths.

Bandwidths represent the number of nearest county-level observations included in the adaptive bi-square kernel. Smaller bandwidths indicate more localized spatial relationships, whereas values approaching 49 indicate broader or nearly global relationships.

Figures 9 and 10 present the spatial distributions of the MGWR coefficients that passed the variable-specific multiple-testing-adjusted significance threshold. The maps indicate that statistically supported associations were generally broad in spatial extent rather than highly localized, which is consistent with the relatively large variable-specific bandwidths. GDP per capita was significantly negative across all counties in 2010 and 2014, but no county remained significant in 2018 or 2023. Local budget revenue was significantly positive in most counties in 2010 and across the entire study area in 2018, whereas its coefficients were not statistically significant in 2014 or 2023. Nighttime light intensity did not pass the adjusted significance threshold in any of the four selected years.

Road-network density exhibited the most persistent positive pattern, becoming significant throughout the study area in 2014, 2018, and 2023. The 2023 road-network-density result should nevertheless be interpreted cautiously because its annual VIF marginally exceeded 10 in that cross-section. Per capita savings was also significantly positive across all counties in 2010, 2014, and 2018, but became insignificant in 2023. Distance to Kunming displayed a similarly broad positive association during the first three selected years, while significant positive coefficients were restricted to four northeastern counties in 2023. These results indicate that infrastructure, financial capacity, and accessibility were associated with CCD-R over relatively broad spatial scales, although the statistical support for savings and distance weakened in the final year.

The cultural variables showed considerably weaker and less stable associations. ICH density was not significant in 2010, 2018, or 2023, but displayed significant negative coefficients in 41 counties in 2014. Minority population share was significantly negative in only two southern counties in 2010 and significantly positive in only one northern county in 2018; no significant association was identified in 2014 or 2023. Overall, the number of statistically supported associations narrowed over time, with road-network density emerging as the only factor that remained significant across all counties in 2023. The findings therefore suggest that spatial heterogeneity was selective rather than pervasive and that cultural characteristics had comparatively weak and temporally unstable local associations with CCD-R.

4.4.3 XGBoost performance and SHAP results

The XGBoost model demonstrated stable predictive performance under the three leakage-controlled validation strategies. In the nested spatial GroupKFold validation, in which counties in the test folds were excluded from model training and hyperparameter selection, the pooled out-of-fold predictions achieved an R² of 0.763, an RMSE of 0.0526, and an MAE of 0.0400. The mean fold-level R² was slightly lower at 0.718, indicating some variation in predictive performance across different groups of held-out counties. Under the expanding-window temporal validation, the model achieved a pooled R² of 0.777, an RMSE of 0.0482, and an MAE of 0.0401 when predicting 2018–2023 using only observations from preceding years. The strict spatiotemporal holdout, which simultaneously excluded ten counties and later-year observations from training, produced an R² of 0.851, an RMSE of 0.0443, and an MAE of 0.0373. Although the latter result indicates that model performance remained robust when both spatial and temporal overlap were removed, it was based on only 40 test observations and is therefore treated as supplementary evidence rather than as a directly comparable improvement over the larger cross-validation samples (Tables 13, 14).

Table 13. Predictive performance of the XGBoost model under alternative validation strategies.

The nested spatial validation assigned counties to non-overlapping training and test folds. The expanding-window validation used only years preceding each test year. The strict spatiotemporal holdout tested later observations from ten counties that were excluded from model training. The strict holdout result is based on a smaller test sample and should be interpreted as supplementary robustness evidence.

Table 14. Out-of-fold SHAP importance of the explanatory variables.

SHAP importance was calculated from pooled out-of-fold predictions generated by the nested county-grouped cross-validation models. It represents contribution to model prediction and should not be given a causal interpretation.

The out-of-fold SHAP results identified local budget revenue and per capita savings as the two most influential predictors of CCD-R, accounting for 28.43% and 23.43% of the total mean absolute SHAP contribution, respectively. Together, these two variables represented 51.86% of the model’s total feature importance. Road-network density ranked third, contributing 13.77%, followed by nighttime light intensity at 12.63%. These four fiscal, financial, and infrastructure-related variables jointly accounted for approximately 78.25% of the overall predictive contribution. The remaining factors had more moderate importance: intangible cultural heritage density accounted for 7.91%, GDP per capita for 6.28%, distance to Kunming for 4.63%, and minority population share for 2.93%. The two cultural variables together contributed 10.84%, indicating that cultural characteristics provided additional predictive information but were less influential than fiscal capacity, household savings, transport infrastructure, and urban activity.

The SHAP distribution further revealed clear differences in the direction and functional form of these associations (Figure 11). Higher local budget revenue and per capita savings were generally associated with positive SHAP values, whereas observations with lower values tended to reduce predicted CCD-R. Road-network density showed a similarly positive pattern: counties with sparse road networks were predominantly associated with negative contributions, while the contribution increased rapidly at intermediate levels before gradually approaching a plateau. Nighttime light intensity also displayed a predominantly positive nonlinear association, with its contribution increasing more strongly after passing from low to moderate levels. GDP per capita shifted from negative to positive contributions as income increased, although its overall contribution was smaller than that of local budget revenue and savings. Intangible cultural heritage density showed limited contributions at very low levels but became increasingly positive once heritage density reached a moderate level. By contrast, minority population share had the lowest overall importance and generally exhibited a weak negative association at higher values. Distance to Kunming displayed a shallow nonlinear pattern, suggesting that its predictive contribution was not adequately represented by a simple linear distance-decay relationship.

The LOWESS-smoothed SHAP curves and county-grouped bootstrap procedure identified several statistically stable observed turning points, as summarized in Table 15 and visualized in Figure 12. Local budget revenue changed from a predominantly negative to a positive predictive contribution at approximately CNY 1.027 billion, with a 95% bootstrap interval of approximately CNY 0.989–1.069 billion. The corresponding turning point for per capita savings was approximately CNY 39,963 per person, with a comparatively narrow interval of CNY 38,777–41,599 per person. GDP per capita crossed the zero-SHAP line at approximately CNY 51,381 per person, while nighttime light intensity and road-network density exhibited turning points of approximately 26.23 and 0.267, respectively. The observed turning point for intangible cultural heritage density was approximately 0.00095. These estimates indicate where the modelled contribution changed direction within the observed sample and should not be interpreted as universal policy thresholds or causal treatment levels. The turning-point intervals for distance to Kunming and minority population share were considerably wider, particularly for distance, and these estimates should therefore be interpreted more cautiously.

Table 15. Observed turning points derived from the SHAP response curves.

Overall, the leakage-controlled XGBoost and out-of-fold SHAP results suggest that variation in CCD-R was most strongly associated with fiscal capacity, household financial resources, transport connectivity, and urban economic activity. Cultural heritage density provided a smaller but non-negligible positive contribution at moderate and higher values, whereas minority population share had only limited independent predictive importance after the other socioeconomic and accessibility variables were considered. These findings complement the panel and MGWR analyses by revealing nonlinear and potentially non-monotonic predictive relationships, while avoiding causal interpretations of the machine-learning results. Supplementary Figures S2–S4 provide additional descriptive spatial context by comparing the distribution of coupling coordination with selected natural conditions, production–living–ecological space dynamics, and cultural characteristics. Because these maps indicate spatial co-occurrence rather than model-based statistical relationships, they are interpreted descriptively only.

4.4.4 Robustness and sensitivity

The panel-regression robustness results are reported in Table 16. The estimates showed some sensitivity to model specification. Road-network density remained positive but statistically insignificant across all alternative models. Per capita savings was negative and significant in the baseline and lagged CCD-R specifications, but its statistical significance weakened under winsorization, alternative standard errors, exclusion of 2020–2022, and alternative dependent-variable specifications. The leave-one-prefecture-out analysis nevertheless showed complete directional consistency, with road-network density remaining positive and per capita savings remaining negative in every subsample. These results support the relative stability of coefficient direction, but not uniform robustness in statistical significance.

Table 16. Panel-regression robustness results.

All models include county and year fixed effects. Baseline standard errors are clustered at the county level. Continuous explanatory variables were standardized, so coefficients represent changes associated with a one-standard-deviation increase.

The MGWR sensitivity tests showed that excluding minority population share, or both cultural variables, generally preserved the direction and spatial distribution of the main local coefficients, although the magnitude and significance of several variables changed in individual years. The leakage-controlled XGBoost results were also stable, with R² values of 0.763 for spatial GroupKFold validation, 0.777 for temporal validation, and 0.851 for the strict spatiotemporal holdout. Overall, the CCD patterns were broadly stable across alternative indicator compositions and subsystem contribution coefficients, while the spatial-autocorrelation results were stable across alternative spatial-weight matrices. Alternative indicator-weighting methods were not evaluated.

5 Discussion

5.1 Spatiotemporal differentiation and the transition in coupling coordination

The results reveal clearly divergent trajectories between socio-economic development and ecological health in the Central Yunnan Urban Agglomeration. Socio-economic development increased rapidly across all five prefecture-level divisions, but this growth did not eliminate the established regional hierarchy. Kunming remained the dominant development core throughout 2010–2023, while Qujing and Honghe gradually emerged as secondary centers. Although peripheral counties experienced substantial relative growth, the absolute gap between the highest- and lowest-performing regions continued to widen. The decrease in the Global Moran’s I of socio-economic development after 2017 suggests that medium-level development gradually diffused beyond the central core; however, the persistence of significant positive spatial autocorrelation indicates that this diffusion remained incomplete. The regional structure therefore evolved from a predominantly low-level and highly concentrated pattern towards a hierarchical system comprising a dominant core, several secondary centers, and improving but still disadvantaged peripheral areas.

Ecological health followed a much less pronounced upward trajectory. The regional mean fluctuated around a relatively stable level, increasing moderately before 2017 but declining again in 2023. Spatially, relatively favorable ecological conditions remained concentrated in parts of the central agglomeration, whereas a persistent low-value zone was located mainly in the eastern counties. The sharp increase in Moran’s I in 2022 and 2023 indicates stronger spatial polarization rather than uniform ecological improvement. In other words, counties with favorable ecological conditions became increasingly concentrated near one another, while the eastern low-value area also became more spatially continuous. This pattern demonstrates that ecological-health disparities were more structurally persistent than the differences in socio-economic development.

These contrasting subsystem trajectories are essential for interpreting the improvement in coupling coordination. The mean CCD-B increased from 0.378 in 2010 to 0.516 in 2023, while the proportion of counties in the severe- and mild-imbalance categories declined from 75.51% to zero. Nevertheless, barely coordinated counties still accounted for 61.22% of the study area in 2023. The urban agglomeration had therefore largely moved beyond explicit imbalance, but most counties had not yet achieved high-level coordination. Because socio-economic development increased much more rapidly than ecological health, the improvement in CCD should not be interpreted as evidence that both subsystems advanced at the same rate. Rather, it indicates a gradual transition from low-level imbalance towards intermediate coordination.

The reconstructed CCD-R confirms that this broad conclusion was not produced solely by overlapping socio-economic indicators. Although CCD-R generated more conservative absolute values and altered the classification of some county-year observations, it retained a strong correlation with CCD-B and increased from 0.374 to 0.496 between 2010 and 2023. The central high-value structure also remained visible. This indicates that the long-term improvement and core–periphery hierarchy were broadly robust, while also showing that discrete coordination classes are more sensitive than continuous values and rankings to indicator composition. The supplementary comparisons of natural conditions and production–living–ecological space dynamics provide additional geographic context for these patterns (Supplementary Figures S2, S3), but they are interpreted as descriptive spatial correspondence rather than evidence of independent effects.

5.2 Multiscale and nonlinear associations of socioeconomic and infrastructure factors

The panel, MGWR, and XGBoost–SHAP results provide complementary evidence because they represent different dimensions of association. The fixed-effects models estimate whether changes within the same county over time are associated with changes in CCD-R. The Mundlak model separates these within-county relationships from persistent differences between counties. MGWR examines whether cross-sectional coefficients vary spatially, whereas XGBoost–SHAP identifies nonlinear predictive contributions across the pooled panel. Differences among the results should therefore be understood in relation to these different estimands rather than regarded as direct contradictions.

This distinction is particularly apparent for road-network density. Annual within-county changes in road density were positively but not significantly associated with CCD-R in both the contemporaneous and lagged fixed-effects models. In contrast, the long-term county mean of road-network density was positively associated with coordination in the Mundlak model. Road density also ranked third in SHAP importance and was the only MGWR variable with significant positive coefficients across all counties in 2023. However, the 2023 local coefficient pattern should be interpreted cautiously because road-network density showed a VIF slightly above 10 in that year. Taken together, these findings suggest that the accumulated availability of transport infrastructure is more closely associated with socio-ecological coordination than short-term annual changes in road construction. Transport infrastructure may generate its main benefits through long-term improvements in accessibility, public-service provision, and regional connectivity rather than through immediate year-to-year changes.

Per capita savings exhibited a similarly important within–between distinction. In the fixed-effects models, within-county increases in savings were negatively associated with CCD-R, and the lagged coefficient remained significant. However, counties with higher long-term average savings had higher coordination levels in the Mundlak model, while savings ranked second in the SHAP importance analysis. The negative within-county coefficient was also sensitive to alternative specifications and should therefore not be interpreted as evidence that household savings reduce coordination. It may instead capture short-term adjustment processes or unobserved county-level shocks that cannot be fully identified in the current design. The contrast demonstrates that temporary changes within a county and persistent structural differences across counties represent substantively different relationships.

Fiscal capacity emerged as the strongest predictive dimension. Local budget revenue accounted for 28.43% of the total mean absolute SHAP contribution and showed comparatively localized MGWR bandwidths in 2010 and 2014. This suggests that fiscal resources were associated not only with general regional differences but also with more localized patterns during the earlier study period. Fiscal capacity may reflect the ability of local governments to support infrastructure, welfare provision, ecological management, and public services, although the present results do not identify which expenditure channel is responsible for the observed association.

Nighttime light intensity and GDP per capita also displayed nonlinear patterns. Their SHAP contributions generally changed from negative or weak values at lower levels to positive values at moderate and higher levels. Road-network density showed a rapid increase in predictive contribution before approaching a plateau. These findings indicate that the relationships between development conditions and CCD-R cannot be fully represented by single linear coefficients. However, the estimated zero crossings should be interpreted as turning points observed within the study sample rather than as universal policy thresholds. Their locations may vary across regions, periods, measurement systems, and model specifications.

The MGWR diagnostics further qualify the interpretation of spatial heterogeneity. MGWR achieved relatively high explanatory power, but it did not consistently outperform OLS or conventional GWR according to AICc. Moreover, most variable-specific bandwidths included 45–47 of the 49 counties, indicating that many relationships operated at broad regional or near-global scales. Stronger local variation was concentrated mainly in local budget revenue during the earlier years and in a limited number of cultural coefficients. MGWR is therefore most appropriately regarded as a supplementary tool for identifying selected local differences rather than evidence that every explanatory relationship varies strongly across space.

5.3 Cultural context and differentiated governance implications

The revised empirical results support a more cautious interpretation of cultural characteristics. ICH density accounted for 7.91% of total SHAP importance and showed a positive predictive contribution at moderate and higher observed levels. Minority population share, by contrast, contributed only 2.93% and showed limited or unstable MGWR significance. The two variables together added predictive information, but their combined contribution remained substantially smaller than that of local budget revenue, savings, road-network density, and nighttime light intensity.

These findings indicate that ethnic culture should be interpreted as a contextual correlate rather than as a dominant explanatory factor of coupling coordination. ICH density measures the spatial concentration of officially recognized heritage resources, while minority population share is a broad demographic measure. Neither variable directly captures cultural vitality, traditional ecological knowledge, local participation, or the effectiveness of heritage governance. The descriptive comparison in Supplementary Figure 4 is therefore retained only to illustrate the cultural geography of the study area. Formal interpretation should rely on the MGWR and out-of-fold SHAP results rather than on visual overlap between cultural maps and high-CCD areas.

The findings nevertheless have important implications for differentiated regional governance. In Kunming’s central districts, where socio-economic and coordination levels were already relatively high, the primary challenge is to improve the ecological quality and efficiency of growth. Policy should focus on controlling land-intensive expansion, improving resource-use efficiency, strengthening ecological restoration, and ensuring that fiscal and infrastructure investment produces ecological and welfare benefits rather than economic output alone.

For secondary urban centers in Qujing, Yuxi, Chuxiong, and Honghe, the results support continued improvement in transport connectivity, fiscal capacity, and public services, but such investment should be coordinated with industrial upgrading and ecological carrying-capacity constraints. The SHAP response patterns suggest that the contribution of several development variables weakens after moderate levels are reached. Expanding infrastructure or economic activity without improving development quality may therefore yield limited additional coordination gains.

Peripheral counties require greater attention to basic accessibility, social welfare, and locally appropriate livelihood opportunities. At the same time, the persistent eastern ecological low-value zone indicates that ecological restoration and environmental management should extend beyond individual county boundaries. Cultural heritage may be incorporated into these strategies through community-based conservation, locally appropriate rural development, and low-impact tourism. However, cultural-resource abundance should not be assumed to improve coordination automatically; its contribution depends on whether cultural assets are translated into household welfare, environmental stewardship, and effective local governance.

5.4 Methodological implications, limitations, and future research

The revised analytical framework offers several methodological improvements. First, pooled entropy weights were held constant across the complete study period, strengthening intertemporal comparability. Second, the reconstructed CCD-R reduced mechanical overlap between the dependent variable and the explanatory variables. Third, the panel models distinguished within-county changes from between-county structural differences. Finally, the XGBoost analysis used spatial GroupKFold, expanding-window temporal validation, and strict spatiotemporal holdout, while SHAP values were calculated from out-of-fold predictions. These modifications allow descriptive patterns, statistical associations, and predictive contributions to be distinguished more clearly without assigning them a causal interpretation.

The robustness results also show that stability is not uniform across all findings. The continuous CCD values, county rankings, and broad spatiotemporal pattern remained stable under alternative index specifications, but coordination classifications were more sensitive to joint indicator exclusion. Similarly, the direction of the panel coefficients was generally stable, whereas their statistical significance changed across model specifications. Excluding cultural variables preserved most MGWR coefficient patterns, although the magnitude or significance of some local estimates changed. Robustness should therefore be understood as broad consistency in the main regional patterns rather than identical results under every specification.

Several limitations remain. The county-level scale conceals substantial variation between urban centers, peri-urban zones, rural settlements, and ecologically sensitive areas. The ecological-health index emphasizes vegetation productivity, forest cover, landscape diversity, and resilience, but provides limited representation of pollution pressure, water quality, climate exposure, and direct ecological degradation. The equal subsystem weights and coordination-class thresholds also involve normative choices, even though they are commonly used in coupling-coordination research. Alternative indicator-weighting methods were not examined in this study; therefore, the robustness results should not be interpreted as demonstrating insensitivity to the entropy-weighting scheme.

The measurement of cultural characteristics remains particularly constrained. Minority population share was derived primarily from census observations and approximated for non-census years, while ICH density does not measure the actual intensity of cultural practice or traditional ecological knowledge. The annual MGWR models also contain only 49 observations, and most selected bandwidths approach the full sample size. Furthermore, the observed SHAP turning points are specific to the present data and predictive model.

Future research should therefore incorporate finer spatial units, richer environmental-monitoring data, and more direct indicators of cultural participation, heritage governance, and traditional ecological practices. And future assessments could also distinguish regional status, progress, and distance from desired sustainability conditions and compare fixed pooled entropy weights with alternative goal-oriented or profile-based weighting approaches (Ma et al., 2026). Alternative weighting methods, subsystem weights, and coordination classifications should also be compared systematically. Longer panels, external validation across other mountainous urban agglomerations, dynamic spatial models, and quasi-experimental designs would help determine whether the observed associations are generalizable across contexts and, where appropriate, support stronger causal inference.

6 Conclusion

This study examined the spatiotemporal evolution and coupling coordination of socio-economic development and ecological health across 49 county-level units in the Central Yunnan Urban Agglomeration from 2010 to 2023. A pooled entropy-weighted TOPSIS framework was combined with spatial autocorrelation analysis, alternative coupling-coordination measures, panel regression, MGWR, and leakage-controlled XGBoost–SHAP analysis. To address the potential circularity between the constructed dependent variable and the explanatory variables, a strict non-overlapping coordination measure, CCD-R, was reconstructed and used as the principal outcome in the associated-factor analyses.

The results yield four main conclusions. First, socio-economic development increased rapidly but retained a persistent core–periphery structure centered on Kunming, while secondary development nodes gradually emerged in Qujing, Honghe, Yuxi, and Chuxiong. Ecological health remained comparatively stable and spatially differentiated, with a central relatively favorable zone and a persistent eastern low-value cluster. Second, coupling coordination improved substantially. The regional mean CCD-B rose from 0.378 to 0.516, and counties classified as severely or mildly imbalanced declined from 75.51% in 2010 to zero in 2023. Nevertheless, barely coordinated counties still constituted the majority in 2023, indicating that regional coordination had improved but had not reached a uniformly high level. The reconstructed CCD-R preserved the broad temporal trend and spatial hierarchy, confirming that the principal findings were not solely generated by overlapping indicators.

Third, the associated factors operated differently across temporal, spatial, and predictive dimensions. Within-county changes in road-network density were not significantly associated with CCD-R, whereas long-term differences in road density were positive. Per capita savings displayed a negative within-county association but a positive between-county and nonlinear predictive relationship. The MGWR results revealed localized variation for selected fiscal and cultural variables, although most relationships operated at broad spatial scales and MGWR did not consistently outperform simpler models. Fourth, the XGBoost–SHAP analysis identified local budget revenue, per capita savings, road-network density, and nighttime light intensity as the most influential predictors. Cultural characteristics contributed additional predictive information, but ICH density had only moderate importance and minority population share had limited independent importance. The evidence therefore supports interpreting cultural characteristics as contextual correlates rather than dominant explanatory factors.

Overall, the study shows that improving socio-ecological coordination in mountainous urban agglomerations requires more than accelerating economic growth. Policy should focus on converting fiscal, financial, and infrastructure capacity into ecological and welfare improvements, strengthening secondary development centers without reproducing excessive environmental pressure, and addressing persistent ecological disadvantages in peripheral areas. Cultural heritage can support locally adapted development, but its contribution depends on broader institutional, economic, and ecological conditions. Because the reported relationships are associational and the observed turning points are sample-specific, the findings should inform differentiated governance rather than be given a causal interpretation or treated as fixed policy thresholds.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2026.1897877/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Purabi Saikia, Banaras Hindu University, India

Reviewed by: Qing Lu, Chinese Research Academy of Environmental Sciences, China

Yujia Chen, Tianjin University, China