## 4 Materials and Methods  
The data and methodology are chosen to enable testing of the SATUM model and correspond to research questions 1-4. For this study, nearly all of the provincial data were drawn from the Chinese Statistical Yearbooks in 2021 and 2016. The Chinese Statistical Yearbook is an annual publication compiled by the National Bureau of Statistics of China and serves as an authoritative source of official statistical information [37]. It has been widely used for research [46; 48; 49]. Although some weaknesses have been identified, namely problems with time series industrial classification definitions over time, coverage of the industrial sector [50], difficulties in searching for terms in English [48], and questions about data collection of certain variables [45], we feel the present dataset has limited exposure to these problems based on the detailed notes in the yearbooks. We do not include GDP or industrial production data, which have been criticized [45; 48; 50]. The variables are either indices or normalized by population, and all are from the years 2020 and 2015, except for unemployment, which was averaged for the year and the two preceding ones in order to even out year-to-year fluctuations. The data were checked for accuracy and entered into SPSS (v26.0) for statistical analysis and ArcGIS Pro (v3.2) for GIS mapping. First, the Moran's I diagnostic for each of the dependent variables was calculated to specify the amount of spatial agglomeration or dis-agglomeration. Moran's I varies from 1.0 (nearby provinces are similar to each other) to \(-1.0\) (nearby provinces are different from each other). A value of zero indicates that the digital divide variables are independently and randomly distributed in space. Moran's I was also calculated for the regression residuals to assess the independence and randomness of residuals. A significant Moran's I value for residuals indicates that the regression findings must be viewed with caution, discarded, or re-run with new variables that account for the bias. Besides descriptive statistics, ordinary least squares (OLS) regression was employed to test proposed associations between the 16 dependent variables and 14 independent variables. OLS regression was chosen instead of alternative approaches, such as structural equation modeling [49; 51; 52; 53] or path analysis [54] for the following reasons. For SEM, there are modeling assumptions that need to be met, including a cause-and-effect relationship such that cause precedes effect in time, linearity, and other assumptions [52]. The main reason is that the present sample size of 31 is too small to apply structural equation modeling. Studies have shown that SEM sample sizes commonly range from 200 to 400 but are certainly more than 50 [49; 52; 55]. Likewise, another type of multi-factor modeling, path analysis, requires sample sizes that are too large for the present sample size of 31. [PERSON] [54] indicates a minimum sample size for path modeling of 10 times the number of parameters. The present study has 4 or more parameters, implying a minimum sample size of 40 or more. OLS regression was conducted in a stepwise manner, allowing variables to enter with significance levels equal to or less than 0.10. The Variance Inflation Factor (VIF) was computed for each independent variable as a test of multicollinearity, with a cutoff of 5.0 [56]. Each of the 14 OLS regressions had a VIF lower than 5.0, indicating that multicollinearity problems were not present. All OLS regression models of this study were computed using IBM SPSS software, version 29, which produced diagnostics, such as the VIF values of each OLS model. OLS regression diagnostics of Joint Wald, Koenker, and Jarque-Bera were calculated to ensure that the assumptions of OLS regression are met (see discussion in [25], (pp. 9-11)). These diagnostics were obtained from the spatial statistics module of the ArcGIS Pro software, which is often used for mapping and geospatial analysis.