Socioeconomic and Environmental Impacts on Regional Tourism across Chinese Cities: A Spatiotemporal Heterogeneous Perspective

Abstract

Understanding geospatial impacts of multi-sourced drivers on the tourism industry is of great significance for formulating tourism development policies tailored to regional-specific needs. To date, no research in China has explored the combined impacts of socioeconomic and environmental drivers on city-level tourism from a spatiotemporal heterogeneous perspective. We collected the total tourism revenue indicator and 30 potential influencing factors from 343 cities across China during 2008–2017. Three mainstream regressions and an emerging local spatiotemporal regression named the Bayesian spatiotemporally varying coefficients (Bayesian STVC) model were constructed to investigate the global-scale stationary and local-scale spatiotemporal nonstationary relationships between city-level tourism and various vital drivers. The Bayesian STVC model achieved the best model performance. Globally, eight socioeconomic and environmental factors, average wage (coefficient: 0.47, 95% credible intervals: 0.43–0.51), employed population (−0.14, −0.17–−0.11), GDP per capita (0.47, 0.42–0.52), population density (0.21, 0.16–0.27), night-time light index (−0.01, −0.08–0.05), slope (0.10, 0.06–0.14), vegetation index (0.66, 0.63–0.70), and road network density (0.34, 0.29–0.38), were identified to have nonlinear effects on tourism. Temporally, the main drivers might have gradually changed from the local macro-economic level, population density, and natural environment conditions to the individual economic level over the last decade. Spatially, city-specific dynamic maps of tourism development and geographically clustered influencing maps for eight drivers were produced. In 2017, China formed four significant city-level tourism industry clusters (hot spots, 90% confidence), the locations of which coincide with China’s top four urban agglomerations. Our local spatiotemporal analysis framework for geographical tourism data is expected to provide insights into adjusting regional measures to local conditions and temporal variations in broader social and natural sciences.

1. Introduction

Tourism is an underlying industry that promotes the development of the global economy [1]. According to the World Travel & Tourism Council (WTTC), tourism contributed 10.3% (8.9 trillion US dollars) of global GDP and provided one-tenth of the total number of jobs (330 million positions) in 2019 before the pandemic [2]. Through developing the tourism industry, local governments can markedly improve the level of infrastructure construction, increase employment opportunities, improve people’s living conditions, and promote urban economic growth [3,4,5]. In addition, tourism development is a fundamental part of a sustainable development strategy, which is recognized as a green industry by the world due to its low energy consumption and light pollution characteristics in the development process [6].

Despite being one essential force promoting regional economy, regional tourism itself is greatly influenced by socioeconomic status [7,8], including GDP [8], employment status [9], personal income [10], health and hygiene [7], industrial production index [11] and social media [12]. Besides the socioeconomic condition, research also identified the notable role of the environment in affecting regional tourism [13,14,15], especially climatic conditions, such as temperature [16], precipitation [17], sunshine [18], and relative humidity [19]. Road infrastructure was also a critical environmental driver enhancing the tourism industry [20,21]. However, all these previous studies only adopted a limited number of factors. It is necessary to consider the comprehensive impacts on tourism by combining socioeconomic conditions with environmental conditions.

When investigating relationships between regional tourism and potential explanatory factors, an unrealistic assumption persistently embedded in previous literature was that the variables’ relationships were homogeneous, which had been defined as stationarity. For instance, non-spatial tourism studies using qualitative analysis [22], feasible generalized least square (FGLS) regression [19], linear and quantile regression [23], or logit regression [24] are regarded as global-scale analyses and also ignore the existence of spatial effects. Likewise, some geospatial tourism studies using the spatial regression models, such as the exploratory spatial data analysis (ESDA) [25] or spatial econometric models [26], are capable of incorporating spatial effects for intercept or residual but are still unable to estimate a set of space-scale coefficients to characterize the varying region-specific relationships between variables. Hence, a more reasonable assumption in the real world highlights the heterogeneous or varying impacts of explanatory variables on tourism development due to region-specific situations, especially for studies conducted across large domains at finer geospatial scales. Such spatially heterogeneous variables relationships are called spatial nonstationarity in the field of statistics. At present, the geographically weighted regression (GWR) [27] is frequently used in tourism research, aiming at exploring such spatial nonstationarity between tourism and various influencing factors [28,29]. However, to the best of our knowledge, no study has been conducted from the spatiotemporal integrated nonstationary perspective, to fully explore both socioeconomic and environmental drivers on regional tourism development.

In China, as the area of interest in this study, there has long been an issue of regional tourism development disparities [30], which obstructed regional tourism sustainability to some extent [31]. Although these geospatial disparities have been extensively discussed at a provincial-level scale [32] or city group scale [33], seldom have studies explored the city-specific disparities of regional tourism, especially over mainland China. Based on tourism connotations and tourism elements, Chinese scholars have established a comprehensive indicator framework of influencing the urban tourism industry from multiple dimensions. Socioeconomic and environmental aspects are also considered indispensable indicators reflective of regional tourism industry development [34]. However, no existing studies ever investigated the joint impacts of socioeconomic and environmental conditions on China’s city-level tourism from a spatiotemporal heterogeneous perspective, to provide evidence-based implications for assisting the formulation of tourism-related policies at governmental levels in a timely and effective manner.

In an attempt to find effective factors affecting regional tourism outcomes to provide tourism strategies tailored for specific local spatial conditions and changing temporal circumstances, we constructed an explanatory variable framework composed of 30 variables, including socioeconomic and environmental conditions. We explored spatiotemporal heterogeneous relationships between the regional tourism development and the multi-source explanatory factors from 2008 to 2017 across Chinese cities by employing the Bayesian spatiotemporally varying coefficients (STVC) model [35,36]. The establishment of such an explanatory variable framework in our study also served as a contributor to the current literature in this field in terms of improving the comprehensiveness of the existing research index system as well as adding novel perspectives into this field based on the consideration of both spatial and temporal heterogeneity.

2. Materials and Methods

2.1. Study Area and Data

Considering the unbalanced development speed and regional differences in China’s tourism industry during the last decade, in this study, 343 prefecture-level areas were selected as the underlying research units (excluding Hong Kong, Macao, and Taiwan). Total tourism revenue was employed as a proxy variable to describe the regional tourism development level from 2008 to 2017 [30]. Figure 1 illustrates the original geographical distribution of city-level total tourism revenue across China in 2017.

Correspondingly, we collected a relatively comprehensive system of 30 explanatory variables at the city level, including 21 socioeconomic factors and nine environmental variables (summarized in

Table 1. City-level potential explanatory variables of regional tourism in China: SV1-21 denote twenty-one socioeconomic factors, and EV1-9 denote nine environmental factors.

Identifier	Socioeconomic Variables	Identifier	Environmental Variables
SV1	Gross domestic product (GDP) per capita (yuan)	EV1	Precipitation (0.1 mm)
SV2	Population density (person/km2)	EV2	Temperature (centigrade)
SV3	Employment density of the first industry (person/km2)	EV3	Air pressure (1 N/m2)
SV4	Employment density of the Second industry (person/km2)	EV4	Humidity (hPa)
SV5	Employment density of the tertiary industry (person/km2)	EV5	NDVI (/)
SV6	Mobile phone penetration rate (subscriber/person)	EV6	Road network density (km/km2)
SV7	Internet broadband penetration rate (subscriber/person)	EV7	Elevation (meter)
SV8	Local general budget revenue per capita (yuan)	EV8	Slope (°)
SV9	Local government budgetary expenditures per capita (yuan)	EV9	Nighttime light index (/)
SV10	Employees population density (person/km2)
SV11	Savings deposits of per capita residents (yuan)
SV12	Loans of financial institutions per capita (yuan)
SV13	Industrial enterprises density (number/km2)
SV14	Social fixed asset investment per capita (yuan)
SV15	Social consumable total retail sales per capita (yuan)
SV16	Student’s density of ordinary middle school (person/km2)
SV17	Student’s density of primary school (person/km2)
SV18	Hospital density (number/km2)
SV19	Hospital beds per capita (number/person)
SV20	Employment density of urban units (person/km2)
SV21	Average wage of employed persons in urban units (yuan)

), to detect their impacts on total tourism revenue in China. The total tourism revenue and socioeconomic data were retrieved from the China City Statistical Yearbook and Statistical Bulletin. The climate data (EV1-EV4) were collected from the National Meteorological Information Center (http://data.cma.cn/, accessed on 28 April 2021). The other environmental factors (EV5-EV8) were downloaded from the Resource and Environment Science and Data Center (http://www.resdc.cn/, accessed on 28 April 2021). As a list of environmental variables, including elevation, road network density, slope, and nighttime light index, were not temporally continuous data, these variables were only added as a part of the local-scaled modeling for spatial nonstationary analysis. Other socioeconomic and environmental factors had spatiotemporal variation characteristics, satisfying the hypothesis of spatiotemporal nonstationarity.

2.2. Statistical Methods

2.2.1. Variable Selection

Two widely adopted approaches, namely multicollinearity assessment and random forest [37], were employed in a progressive manner as a screening step for identifying the most representative influencing factors on the tourism industry from 30 candidate variables. Precisely, the indicator variance inflation factor (VIF) was first adopted to measure the multicollinearity effect, referring to a correlation between explanatory factors [38]. Commonly, VIF < 10, representing mild and negligible multicollinearity, is adopted as the threshold to screen variables [39]. Here, given the adequacy of candidate variables involved in this analysis, a stricter standard was adopted, indicating that a candidate variable with VIF > 5 was removed. Following the VIF step, random forest, an integrated machine learning approach relying on the decision tree, was adopted for further screening the explanatory variables according to the calculation of an indicator named mean decrease impurity (MDI), which has been commonly used for reflecting the ranking of a factor’s relative importance [40]. For a candidate variable, a higher value of MDI is associated with the increased importance of the variable. This random forest step is typically empirical and data-driven, as MDI is not a relative measure [41].

2.2.2. Bayesian STVC Model

The Bayesian spatiotemporally varying coefficients (STVC) model is a recently burgeoning local spatiotemporal regression developed under the Bayesian hierarchical modeling (BHM) framework. It is mainly designed to quantitatively characterize structured and heterogeneous spatiotemporal impacts (expressed as local-scale coefficients) of different covariates on the outcomes of the variable of interest, that is, to explore the spatiotemporal nonstationarity inherent in geospatial research phenomena [35,36].

For China’s tourism case, $Y_{it}$ denotes the space–time monitoring data of the total tourism revenue indicator, in which $i=1,\dots ,I$ (I = 343) are the administrative geographical units of the cities. For each city, data are available for a ten-year period from 2008 to 2017, labeled as $t=1,\dots ,T$ (T = 10). Then, the structured additive predictor ${\zeta }_{it}=g(Y_{it})$ within a reduced Bayesian STVC model is formulated in Equations (1)–(3), i.e.,

$${\zeta }_{it}=g(Y_{it})=\eta +{\displaystyle {\sum }_{k=1}^K{f}_{space}({\omega }_{ik}^{}S{X}_{itk})}+{\displaystyle {\sum }_{m=1}^M{f}_{time}({\phi }_{tm}^{}T{X}_{itm})},$$

(1)

$${\omega }_i|{\omega }_{-i},{\tau }_{\omega }, W~N(\frac{{\displaystyle {\sum }_{j=1}^I{w}_i{}_j{\omega }_j}}{{\displaystyle {\sum }_{i=1}^I{w}_i{}_j}},\frac{1}{{\tau }_{\omega }{\displaystyle {\sum }_{i=1}^I{w}_i{}_j}}),i=1,\dots ,I,$$

(2)

$${\phi }_{t+1}-{\phi }_t|{\tau }_{\phi }~N(0,\frac{1}{{\tau }_{\phi }}),t=1,\dots ,T-1, \mathrm{or}$$

$${\phi }_t-2{\phi }_{t+1}+{\phi }_{t+2}|{\tau }_{\phi }~N(0,\frac{1}{{\tau }_{\phi }}),t=1,\dots ,T-2,$$

(3)

In Equation (1), $g(\cdot )$ denotes a log-Gaussian likelihood function for this case to link $Y_{it}$ and ${\zeta }_{it}$. $\eta$ denotes the intercept with fixed effect. $SX$ signifies $K$ main covariates with the spatial nonstationary assumption. $TX$ represents $M$ main covariates that are assumed to be temporally nonstationary. The parameter ${\omega }_{ik}$ is named as space-coefficients (SCs) and ${\phi }_{tm}$ is named time-coefficients (TCs), which are two fundamental outputs of the STVC model. $f_{space}(\cdot )$ and $f_{time}(\cdot )$ signify the spatial and temporal latent Gaussian models (LGMs) that are used for fitting the random effects of spatial and temporal nonstationarity to estimate local parameters SCs and TCs [42,43].

In Equation (2), on account of the spatial LGM $f_{space}(\cdot )$, the prior intrinsic conditional autoregressive (iCAR) model is adopted for fitting the spatial autocorrelation characteristics that are also called the spatial structured random effects within a BHM [44], where ${\omega }_{-i}$ denotes every spatial unit in $\omega$ apart from the i-th spatial unit, $W=(w_{ij})$ represents the spatial relation matrix in which $w_{ij}=1$ if spatial units i and j are neighbors, e.g., spatial adjacency relations here, and $w_{ij}=0$ otherwise, as well as ${\tau }_{\omega }$ further indicates the precision parameter [45].

In Equation (3), the prior random walk (RW) model is used as the temporal LGM $f_{time}(\cdot )$ to estimate the temporal autocorrelation characteristics of TCs, where the structured temporal random effect of covariates $\phi$ can be a random walk of order one or two, with ${\tau }_{\phi }$ being the precision parameter [46]. The prior RW model of order two is more suitable for the research object with a clear linear time trend, compared with the prior RW model of order one.

2.2.3. Model Implementation and Comparison

To explore both the global homogeneous and local heterogeneous impacts of socioeconomic and environmental factors on city-specific outcomes of China’s tourism, we implemented four types of Bayesian regressions, the multiple linear regression (MLR, model 1), the ordinary generalized additive model (GAM, model 2), the global spatiotemporal regression (model 3), and the Bayesian STVC model (model 4), which belongs to the local spatiotemporal regression family. We chose these models based on the following considerations. First, model 1 and model 2 were traditional mainstream models. We used them to fit the overall linear and nonlinear impacts of covariates on tourism [47]. Then, we used a widely applied spatiotemporal regression (model 3), which mainly served as a spatiotemporal descriptive tool, to depict the original smoothed spatial variations and temporal trends of tourism in China [42]. However, models 1–3 are regarded as the global-based type of regression, meaning that covariate impacts (coefficients) were homogeneous across space and over time [35]. Given this underlying limitation of stationarity, model 4 was finally employed to explore the structured heterogeneous (varying) impacts of covariates at both space and time scales [36].

To be specific, the equation of an MLR (model 1) is given by

$${\zeta }_{it}=g(Y_{it})=\eta +{\displaystyle {\sum }_{k=1}^K{\chi }_k^{}{X}_{itk}},$$

(4)

where ${\chi }_k$ denotes the overall coefficient of the k-th covariate $X$, which qualifies the linear numerical impacts of explanatory factors on $Y_{it}$.

An ordinary GAM (model 2) is formulated as [48]

$${\zeta }_{it}=g(Y_{it})=\eta +{\displaystyle {\sum }_{k=1}^K{f}_{GAM}({\delta }_{hk}{X}_{itk})},$$

(5)

where $f_{GAM}(\cdot )$ denotes the nonparametric smooth function for fitting a set of coefficients ${\delta }_{hk}$ with $h$ groups, representing the numerical nonlinear impacts of the k-th covariate. Unlike model 1, model 2 is useful in identifying response–covariate numerical nonlinear relationships. However, both model 1 and model 2 cannot consider the spatiotemporal effects essential for geospatial analysis.

A global spatiotemporal regression (model 3) can be modeled with [42,46,49]

$${\zeta }_{it}=g(Y_{it})=\eta +{\displaystyle {\sum }_{k=1}^K{\chi }_k{X}_{itk}}+f_{space}({\mu }_i)+f_{time}({\lambda }_t),$$

(6)

where ${\mu }_i$ signifies the space-intercepts (SIs) representing the structured spatial distribution of $Y_{it}$, ${\lambda }_t$ signifies the time-intercepts (TIs) representing the structured temporal trend of $Y_{it}$, LGMs $f_{space}(\cdot )$ and $f_{time}(\cdot )$ are the same as in Equation (1).

Model 4, as fully introduced in Equations (1)–(3), has been as a reduced Bayesian STVC regression by removing the spatiotemporal random effects of intercepts to ensure noticeable variations of both spatial and temporal nonstationary impacts of different explanatory factors on the target response variable [36,50].

Finally, the optimal model from the above four regressions with the best model fitness and predictability was further utilized to estimate the complete spatiotemporal maps of $Y_{it}$.

2.2.4. Model Inference and Evaluation

Alternative regression models were established using the Bayesian statistics based on the advanced hierarchical modeling strategy, that is, a BHM framework. Non-informative priors were assigned to parameters within the BHM to embody the idea of data-driven modeling [47]. The integrated nested Laplace approximation (INLA) algorithm, an approximate Bayesian inference technique, was adopted to estimate these regression models using the R-INLA package in the R environment [51] due to its advantage of producing reliable estimated results with a relatively short computation time [52]. The performances of these alternative regressions are evaluated in terms of three aspects, including the degree of model fitting, model complexity, and predictive ability [46]. Specifically, the deviation information criterion (DIC) [53] and the Watanabe–Akaike information criterion (WAIC) are used for reflecting the degree of fitting of the Bayesian regression, for which a smaller value indicates a better model fit. Likewise, the complexity of the Bayesian regression is evaluated with two indices ($P_{DIC}$ and $P_{WAIC}$) that can be simultaneously obtained via the adoption of both the DIC and WAIC methods, for which smaller values are also reflective of better model performances. In terms of the model predictive power, a logarithmic score (LS) retrieved from the conditional predictive ordinates under a leave-one-out cross-validation is used, with smaller values associated with better predictive capacities [54].

3. Results

3.1. Selected Drivers for Modeling

As indicated in Figure 2, through setting VIF < 5 as the inclusion criteria, potential explanatory variables with higher rankings of MDI were selected from the screening outcomes and were added into the regression modeling. To be specific, first, we removed factors with higher multicollinearity based on the exclusion threshold of 5 for VIF, as shown in Figure 2a. This step left 13 socioeconomic factors (i.e., SV1, SV2, SV3, SV6, SV7, SV9, SV10, SV14, SV17, SV18, SV19, SV20, and SV21) and four environmental factors (i.e., EV5, EV6, EV8, and EV9). Further, using Figure 2b, we selected the top eight factors (i.e., EV9, SV21, EV8, EV5, EV6, SV20, SV1, and SV2), which had relative higher importance (contribution) to the response variable. Because the selection of MDI is generally empirical, here, the main reason for our choice of MDI is that there was an apparent bluff trend between the two factors of SV2 and SV7. The top eight factors covered four socioeconomic factors and four environmental factors, which was ideal for exploring the combined impacts of the above two critical aspects on tourism development. Hence, based on the perspectives above, the screening threshold applicable to this case is MDI > 200. Summing up the above, a core variables system particularly applicable to China’s tourism case was created, which contained a total of eight critical factors (renamed as X1–X8 in Table 3), and was further incorporated into the next-step regression analysis.

3.2. Model Assessment and Comparison

We assessed the four types of comparative Bayesian regression models’ performances by jointly considering model fitness, complexity, and predictive power, for which a total of five representative evaluation indicators are summarized in

Table 2. Bayesian modeling evaluations of the alternative regressions for China’s tourism case account for model fitness, complexity, and predictive power.

Index	DIC	${\mathit{P}}_{\mathit{D}\mathit{I}\mathit{C}}$	LS	WAIC	${\mathit{P}}_{\mathit{W}\mathit{A}\mathit{I}\mathit{C}}$
Model 1	40,238.39	9.02	5.97	40,262.83	30.05
Model 2	39,153.72	64.89	5.81	39,157.01	66.83
Model 3	31,825.74	348.95	4.72	31,848.15	338.66
Model 4	30,307.19	899.21	4.52	30,345.55	761.76

Model 1–4: multiple linear regression, generalized additive model, global spatiotemporal regression, and local spatiotemporal regression STVC model; DIC: deviance information criterion; WAIC: Watanabe–Akaike information criterion; $P_{DIC}$: effective number of parameters from DIC; $P_{WAIC}$: effective number of parameters from WAIC; LS: logarithmic score.

. Model 4 (STVC) showed the best performance with the minimum assessment indicators DIC, WAIC, and LS. However, for $P_{DIC}$ and $P_{WAIC}$, model 4 demonstrated a notable deficiency and it presented a much higher complexity than the other three mainstream benchmark models. The complexity ($P_{DIC}$) of model 4 was found to be about 99 times higher than that of a multiple linear regression (model 1), which was 2.6 times higher than that of a global spatiotemporal regression (model 3). Two possible reasons were considered for explaining the increased complexity of the STVC model. Specifically, the STVC model demonstrated both superior model fitness and predictive capacity compared with all the other regressions. Moreover, it should be pointed out that the STVC model was the only one that had the capacity for synchronously detecting both temporal and spatial heterogeneous associations between variables to be further interpreted at a space-time scale. Therefore, model 4 (STVC) was selected as the final regression to explore the spatiotemporal heterogeneous relationships between tourism and eight selected explanatory variables, which was also used for producing a series of estimated spatiotemporal distribution maps reflective of the city-level tourism revenue in China.

3.3. Global-Scale Impacts of Drivers

Two kinds of overall impacts of socioeconomic and environmental variables on tourism were estimated: one was the global-scale linear numerical effects based on model 1; the other one was the global-scale nonlinear numerical effects obtained from model 2. We summarized the critical parameters of model 1 in

Table 3. Linear numerical impacts of main drivers on China’s city-level tourism industry.

Variables	Socioeconomic and Environmental Aspects	Mean	SD	Q 0.025	Q 0.975
X1	Average wage of employed persons in urban units	0.4694	0.0212	0.4278	0.5110
X2	Employment density of urban units	−0.1428	0.0170	−0.1761	−0.1095
X3	GDP per capita	0.4660	0.0258	0.4154	0.5165
X4	Population density	0.2136	0.0282	0.1582	0.2689
X5	Nighttime light index	−0.0141	0.0310	−0.0750	0.0467
X6	Slope	0.1013	0.0202	0.0615	0.1410
X7	Normalized Difference Vegetation Index (NDVI)	0.6630	0.0187	0.6263	0.6996
X8	Road network density	0.3382	0.0226	0.2937	0.3826

, including the overall coefficients representing the stationary relationship among variables, standard deviation (SD), and the 2.5% and 97.5% quantiles of Bayesian credible intervals (CIs). In terms of the four socioeconomic variables, X1 and X2 reflected the income level of individual residents, X3 represented the regional macroeconomic development conditions, and X4 represented the population condition. For the other four environmental variables, X5 and X7 represented the city-specific urbanization process and vegetation coverage based on satellite remote sensing data, respectively. X6 and X8 reflected the general geographical situations characterized by topography and transportation, respectively. Except for X2 and X5, the overall coefficients of the other six factors were found to be greater than zero. This finding indicated that most core variables served as positive stimulants for tourism development from a global-scale perspective. Notably, the NDVI (X7), the average wage of employees in urban units (X1), GDP per capita (X3), and road network density (X8) demonstrated more significant impacts on tourism among the eight factors.

Furthermore, the exponent-scale nonlinear numerical effects of the eight selected drivers were illustrated in Figure 3. We noticed that all the variables’ numerical influencing curves had a similar upward trend. At the same time, we identified the varying impacts of each variable across their development process. It is worth mentioning that X2 and X5 were negatively linearly correlated with the tourism industry, which could not be explained directly. While by further analyzing the nonlinear results, only X2 and X5 appeared to have a significant downward trend, leading to the overall negative linear association in Table 3. This finding proved that model 2 had a superior interpretation capacity over model 1 in fitting global-scale numerical impacts. However, both linear and nonlinear numerical modeling results were produced based on a stationary assumption. As a result, these global-scale outputs might smooth or hide the local-scale heterogeneous impacts of different variables on the tourism industry over the entire study area and time frame, particularly for a fine-scale space–time dataset.

3.4. Temporally Varying Impacts of Drivers

In Figure 4, we presented a TIs graph and five TCs graphs with 95% Bayesian CIs, to exhibit the crude temporal dynamic trend of tourism and the temporally heterogeneous impacts of main drivers on tourism in China, as well as the uncertainties of these estimated parameters. According to Figure 4a, China’s tourism development level demonstrated a continuously increasing trend from 2008 to 2017, meaning that China’s tourism industry maintained a high development speed spanning ten years. Furthermore, it can be seen from Figure 4b that the temporal tourism–covariates relationships varied non-linearly over 2008–2017. This visualization of local-scale nonstationary regression relationships over periods was an essential feature of the Bayesian STVC model that cannot be facilitated via the adoption of global-scaled coefficients. Generally speaking, X3 (GDP per capita), X4 (population density), and X7 (NDVI) showed a downward trend from 2008 to 2017, which indicated a strong to weak impact of these variables on tourism development over time. While X1 (average wage of employed persons in urban units) and X2 (employment density of urban units) presented an initial downtrend, followed by an upward tendency starting from 2013, suggesting that their roles in promoting tourism development gradually especially after 2013. These findings also meant that groups with high quality of living might be a potentially vital force to promote tourism. In addition, we noticed that the TCs of X1 and X2 had relatively higher uncertainties (CIs) due to the fluctuation ranges of TCs of X1 (−0.05–0.05) and X2 (−0.02–0.02) being much narrower compared with those of the other indicators, e.g., X4 (−0.4–0.2). In fact, when we plotted all the factors within a single graph using the same vertical coordinate, the uncertainties (CIs) of X1 and X2 turned out to be very small; however, the time trends of X1 and X2 could be smoothed out. From this perspective, the uncertainties of all indicators were within an acceptable range.

3.5. Spatially Varying Impacts of Drivers

Spatially, we retrieved the parameters of SIs from model 3 to geographically map the ten-year average tourism revenue distribution across China, as presented in Appendix A Figure A1. In addition, utilizing the SCs parameters from model 4, the variables’ spatial nonstationary maps were depicted in Figure 5a. The cluster maps for parameter SCs were also produced to highlight those significant (>90% confidence) hot spots and cold spots at the city level, as shown in Figure 5b.

From Figure A1, the spatial distribution characteristics of China’s tourism revenue demonstrated a gradual increase from West China to East China. Simultaneously, we detected diverse geospatial tourism–covariates relationships at the city level from Figure 5a and an apparent spatial agglomeration effect of SC maps in Figure 5b. In fact, for every single factor of interest, city-specific areas with higher sensitivity to this particular covariate could be visually identified in terms of achieving regional tourism development, based on direct analysis of the corresponding SC map produced by that covariate. Furthermore, within each city area, a series of urban policies could be proposed to facilitate its tourism development based on the relative impacts of the eight-core variables. The relative effect of each covariate within each city could be assessed by vertically integrating the local-scale information from all the SC maps together [50].

Looking at the macroscopic regional scale using the hot spot maps in Figure 5b, we may discover that: in Northeast China, X8 may serve as an essential factor for promoting the development of the local tourism industry, while X5 and X6 may have no impacts, and the other factors may also have an individual city-specific impact yet without generating geographic hotspot regions in the past. Likewise, the high-level tourism development in China’s eastern region may be mainly promoted by X2 and X5; and X4 are not entirely essential. In Western China, with low-level tourism development, X4 may be a primary determinant to improve its tourism conditions. Meanwhile, X6 and X7 also present spatially positive clustered effects in some areas, such as Yunnan and Sichuan. In the regions of Central China, tourism development seems to be dominated by socioeconomic factors, including X1, X2, and X3.

3.6. Spatiotemporal Estimated Maps of China’s City-Level Tourism Revenue

A complete series of spatiotemporal distribution maps of China’s city-specific tourism development level from 2008 to 2017 was produced by adopting the optimal Bayesian STVC model (model 4), as shown in Figure 6. The newly model-estimated tourism maps highlighted hidden areas (e.g., cities with missing values) and provided more intuitive information (e.g., smooth the city-level extreme outliers), which were expected to assist in making policies about the sustainable development of tourism. Generally, the overall growth in the city-level tourism industry was identified over the past decade in China, during which time diverse improvement intensities were found among regions at a local city scale. In Central China, since 2008, about 77% of blue-colored cities with weak tourism industries gradually shifted to yellow/red colors with relatively strong tourism industries. In contrast, such a shifting proportion of East China was about 55%, suggesting that Central China’s tourism industry grew faster than that of eastern cities. In terms of West China and Northeast China, the shifting proportions were about 51% and 41%, respectively, which were relatively lower than the other two divisions. In 2017, the low-tourism-level cities with a blue color were mainly distributed in the provinces of Heilongjiang, Gansu, Ningxia, Xinjiang, Qinghai, and Tibet.

Lastly, we performed a hot spot analysis for the newly estimated complete tourism maps in 2007 and 2018, respectively, to detect those significant city clusters (>90% confidence) of the tourism industry, as shown in Figure 7. In 2017, we found four significant tourism industry clusters (hot spots) and one less-developed tourism region (cold spot) at the city level, compared with 2008 with two clearly formed hot spots of the tourism industry. These four identified high-level tourism city clusters in 2017 were demonstrated to be consistent with China’s top four major urban agglomerations, namely, Beijing–Tianjin–Hebei, the Yangtze River Delta, the Pearl River Delta, and the Sichuan–Chongqing Region. This might reveal that current tourism agglomeration development is closely related to the urbanization degree. Meanwhile, a cold spot was detected in West China, indicating that the tourism development of western cities was relatively slow.

4. Discussion

In this study, the multidimensional impacts of socioeconomic and environmental variables, including linear and nonlinear numerical effects and spatiotemporal heterogeneous effects, on regional tourism were comprehensively investigated across Chinese cities along with the first production of a set of spatiotemporal maps depicting China’s total tourism revenue. These findings may add innovative insights about the mechanisms of how multi-source geospatial factors have affected the regional tourism industry, and is expected to provide a brand-new viewpoint for policymakers. According to different scales, we have some conclusions, as follows.

Globally, significant effects of both socioeconomic and environmental variables were identified [28,55,56,57], which highlighted the necessity of taking a wide range of factors into accounts throughout the procedure of tourism policies formulation. Tourism is a comprehensive industry composed of multiple elements, including food, shelter, transportation, travel, entertainment, and purchase. However, the importance of some of these elements embedded in the tourism industry, such as food, shelter, and transportation, is always ignored for the reason that they are simply regarded as the basic service facilities of a city. Therefore, the positive effects of the socioeconomic and environmental factors on tourism are supposed to be focused on the industrial level, which suggests that the idea of developing industries should always be adopted as the guideline for developing the tourism industry regardless of regional or national levels. At present, the “Travel +” strategy being implemented by the Chinese government is exactly based on this idea [58].

Temporally, the development of China’s tourism has mainly benefited from comprehensive time–scale impacts of multiple factors. Based on temporal nonstationarity, the predominant stimulants for tourism development were demonstrated to have gradually switched from the regional economy, populational size, and tourism resource attractiveness to personal economic status. These results implied that China’s current tourism industry demonstrated a new feature that a transition from sightseeing tourism to leisure and holiday tourism is very much likely to occur. Meanwhile, residents’ affluence has been highlighted as an indispensable contributor to nationwide tourism development [59]. Under such a changing background of the tourism industry in China, it is highly suggested that improving personal income, as well as safeguarding the rights and interests of employees, should be adopted as an essential strategy for facilitating the nationwide tourism industry development, which might be achieved via the implementation of multiple tourism-related policies at governmental levels, such as approving paid-leave policies for employees, encouraging enhanced flexibilities of work schedules to be tailored for vocational leaves, as well as encouraging off-peak vocational arrangements.

Spatially, the development of China’s tourism could be characterized as “strong in the east and weak in the west” [30], which was affected by various factors. Cities of West China were mainly affected by population size and tourism resources, while personal income, employment and urbanization had more contributions to cities in the east region [60]. The city-level spatial nonstationarity found in this study could serve as an acceptable reference in the procedure of making more targeted policies by governments at all levels. For example, the western region may put forward corresponding talent introduction policies while promoting economic development. In addition, the local government can develop sightseeing and holiday tourism through developing natural landscapes. Cities of East China need to focus on optimizing the protection system of workers’ rights and interests and developing characteristic tourism products to provide tourists with high-end, comfortable, and personalized services for stimulating tourism. Northeast China may focus on infrastructure and strengthen the planning and laying of the road networks to enhance regional tourism accessibility. Furthermore, city-level local authorities could utilize local resources rationally and determine the direction of tourism strategies by using the critical drivers’ local spatial influencing maps to support ecotourism, sightseeing tourism, vacation tourism, geological tourism, and urban tourism. In addition, the first series of maps displaying China’s tourism revenue’s spatiotemporal distributions at an administrative city level from 2008 to 2017 was produced, which was further analyzed to provide urbanization-related insights into empirically optimizing the unbalanced development of the tourism industry [61,62].

To sum up, from the multidimensional spatiotemporal heterogeneous perspective, the government should formulate various tourism policies based on region-specific conditions, as well as pursue the development concept of “applying proper measurements in line with local conditions and temporal variations”. At present, tourism industry development in areas with relatively high urbanization levels has demonstrated a change from sightseeing tourism to leisure tourism. As a result, socioeconomic status should be continuously considered as a significant factor throughout tourism-related policy-making procedures in these regions. In contrast, regarding cities with low-level urbanization distributed in West China, environmental factors or sightseeing resources, instead of other factors, should be addressed as predominant issues to be considered throughout the formulation of tourism-related policies [60]. Therefore, making city-specific strategies that take city-specific factors into account is expected to improve the accuracy of policy formulation, as well as the effectiveness of strategic implementation, which would further mitigate both “invalid policy” and “weak policy” produced by the “one-size-fits-all” policy.

Finally, we would like to underline the importance of the local spatiotemporal regression approach, namely, the Bayesian STVC model we have selected. As discussed above, introducing a spatiotemporal heterogeneous perspective to regional tourism management could avoid the one-size-fits-all issue via providing multidimensional spatiotemporal information. In the spatial statistics field, local regressions that can deal with such spatiotemporal heterogeneity among variables relationships (spatiotemporal nonstationarity) are relatively rare, which can be generally classified into the frequentist-type model [63,64,65] and the Bayesian-type model [35,36,66], as they were proposed independently under different statistical traditions. The main reasons we chose the Bayesian STVC model as the applied local spatiotemporal regression lie in the following considerations. First, only the Bayesian-based local spatial or spatiotemporal model is an actual “full-map” modeling technique; thus, the results are more reliable [67,68]. Second, the Bayesian STVC model follows a space–time independent nonstationary assumption, dramatically reducing the computational burden and weakening the overfitting problem. Last but not least, due to its separately fitting of space-coefficients (SCs) and time-coefficients (TCs), the Bayesian STVC model is more user-friendly: stakeholders can directly separately obtain the spatial and temporal autocorrelated regularities [36,50]. Beyond these benefits, the Bayesian STVC model still needs further improvement to solve more complex space–time interaction issues in natural and social sciences.

5. Conclusions

This study verifies that socioeconomic and environmental factors simultaneously affect tourism development over China, globally and locally, supported by the up-to-date space-time data of city-level tourism statistics and a series of advanced Bayesian regressions. Remarkably, the local impacts of socioeconomic and environmental conditions vary heterogeneously at the city level in both time and space dimensions across China, and was demonstrated by adopting the cutting-edge Bayesian STVC model, which was also used for estimating the first series of spatiotemporal maps of city-level tourism development. These fruitful findings provide novel insights into policy-making procedures at multiple levels. Here, the Bayesian STVC model was successfully applied to mine the spatial and temporal autocorrelated nonstationarity inherent in tourism–covariates relationships over China and could serve as an emerging tool to offer new insights on spatiotemporal-oriented influencing factor analysis and high-precision prediction in broader GIScience-related fields of social and natural sciences.

Apart from all these achievements, several concerns should be better addressed in future lines of research. First, the seasonal effect is the main factor affecting tourists’ behavior [69], which emphasizes collecting and using quarterly tourism data in tourism research. However, this study is limited because national urban tourism data sources only have annual scale records. Second, other underlying tourism-related factors such as tourism resources were not fully considered in this study [34]. Future studies might focus on a relatively small area with seasonal heterogeneity by using multi-source tourism data to construct more scientific indicators [70] and developing more sophisticated spatiotemporal statistical models for outputting more informative results for regional tourism research.