### Methods

We first construct a healthcare resources supply index (HRS) using the entropy weight method, analyze regional differences, and then examine the sources of differences by using the Gini coefficient decomposition method given by Dagum [36]. Let *G* represent Gini coefficient, *j* and *h* represent the number of regions, *i* and *r* represent the number of cities in the region, *k* represent the total number of regions, *n* represent the total number of cities, and *n*_{j}(*n*_{h}*)* represent the number of cities in region *j*(*h)*. Further, *y*_{ij}(*y*_{rh)} represents the HRS of city *i*(*r)* in region *j*(*h)*, and (overline{y}) represents the arithmetic mean of the HRS. The following equation can be obtained:

$$G{ = }frac{{1}}{{2n^{2} overline{y}}}(sumlimits_{i=1}^{n} {sumlimits_{r=1}^{n} {|y_{i} – y_{r} |} } ) = sumlimits_{j=1}^{k} {sumlimits_{h=1}^{k} {sumlimits_{i=1}^{{n_{j} }} {sumlimits_{r=1}^{{n_{h} }} {|y_{ij} – y_{rh} |/2n^{2} overline{y}} } } }$$

(1)

After the overall Gini coefficient is calculated, *k* regions are sorted according to the average HRS in the YRD, namely *y*_{1} ≤ … ≤ *y*_{j} ≤ …*y*_{k}, and then the Gini coefficient *G* is decomposed into three parts: (1) intra-regional difference contribution (*G*_{w}), (2) inter-regional difference contribution (*G*_{nb}), and (3) transvariation density contribution (*G*_{t}), which meet the requirements of *G* = *G*_{w} + *G*_{nb} + *G*_{t}.

Next, we use the KDE method to explore the spatial distribution dynamics of the HRS in the YRD. Specifically, we use Moran’s *I* to examine the spatial correlation of the YRD’s HRS. If Moran’s *I* is significant, the LISA cluster map will be used to determine whether the local correlation types and clusters in different regions are statistically significant. When Moran’s *I* value is [− 1, 0), 0, and (0, 1], it means negative correlation, irrelevance, and positive correlation, respectively. According to the LISA cluster map, four types of spatial agglomeration relationships can be identified: (1) high-high agglomeration area (HH type), (2) high-low agglomeration area (HL type), (3) low-low agglomeration area (LL type), and (4) low–high agglomeration area (LH type).

Both *σ* convergence estimation and *β* convergence estimation are commonly used in convergence analysis [58]. The former estimate does not rely on an econometric model and uses the coefficient of variation to measure, which can confirm the previous Gini coefficient results. *β* convergence can be divided into absolute *β* convergence and conditional *β* convergence. The former is used to assess whether the supply of healthcare resources in various regions will converge to the same steady-state equilibrium point without considering the city heterogeneity factors. More importantly, we include spatial factors, heterogeneity, and healthcare policy as control variables in conditional *β* convergence for a more comprehensive estimation.

The model for absolute *β* convergence is set as follows:

$$ln (frac{{y_{i,t+1} }}{{y_{i,t} }}) = alpha + beta ln (y_{i,t} ) + mu_{i} + eta_{t} + varepsilon_{it}$$

(2)

where *i* represents the city; *t* represents the time; *y*_{i,t+1} and *y*_{i,t} represent the HRS of city *i* in *t* + 1 period and *t* period, respectively; ln(*y*_{i,t+1}*/y*_{i,t}) represents the annual growth rate of HRS in city *i* during the period from *t* to *t* + *1*; and *β* is the convergence parameter to be estimated, where *β* < 0 means an absolute *β* convergence trend, and otherwise it indicates that there is a divergence trend. *α* is a constant term; *μ*_{i} and *η*_{t} represent regional and time effects, respectively; and *ε*_{it} represents the random interference term. The formula of convergence rate is (v = – frac{1}{TS}ln (1 + beta )), where *TS* represents the time span.

If heterogeneous factors such as economic development status, city size, and government fiscal capacity are included in the model, the model for conditional *β* convergence is set as follows:

$$ln (frac{{y_{i,t+1} }}{{y_{i,t} }}) = alpha + beta ln (y_{i,t} ) + gamma ln (X_{i,t} ) + mu_{i} + eta_{t} + varepsilon_{it}$$

(3)

where *X*_{i,t} is the control variable; *γ* is the parameter to be estimated for the control variable; and the meaning of the other variables is the same as the formula (2).

Considering the possible impact of policies and external shocks on the conditional *β* convergence, we have:

$$ln (frac{{y_{i,t+1} }}{{y_{i,t} }}) = alpha + beta ln (y_{i,t} ) + gamma ln (X_{i,t} ) + delta P_{i,t} + tau [ln (y_{i,t} ) times P_{i,t} ] + mu_{i} + eta_{t} + varepsilon_{it}$$

(4)

where *P*_{i,t} is a policy dummy variable, where if region *i* starts to implement the policy in period *t*, then *P*_{i,t} is assigned a value of 1 during period *t* and after, and *P*_{i,t} is assigned a value of 0 before period *t*; and *δ* represents the impact of the policy on the growth rate of HRS; *τ* represents the impact of policy on the convergence of HRS.

Taking into account the spatial correlation, we construct a spatial econometric model of *β* convergence. In order to determine the exact form of the spatial effects entering the panel model, we first estimate the following three spatial models in the absolute *β* convergence estimation:

$$ln (frac{{y_{i,t+1} }}{{y_{i,t} }}) = alpha + beta ln (y_{i,t} ) + rho sumlimits_{j=1}^{N} {W_{ij} ln (frac{{y_{i,t+1} }}{{y_{i,t} }})} + mu_{i} + eta_{t} + varepsilon_{it}$$

(5)

$$ln (frac{{y_{i,t+1} }}{{y_{i,t} }}) = alpha + beta ln (y_{i,t} ) + mu_{i} + eta_{t} + varepsilon_{i,t} ,varepsilon_{i,t} = lambda sumlimits_{j=1}^{N} {W_{ij} varepsilon_{j,t} } + sigma_{i,t}$$

(6)

$$ln (frac{{y_{i,t+1} }}{{y_{i,t} }}) = alpha + beta ln (y_{i,t} ) + rho sumlimits_{j=1}^{N} {W_{ij} ln (frac{{y_{i,t+1} }}{{y_{i,t} }})} + theta sumlimits_{j=1}^{N} {W_{ij} ln (y_{i,t} )} + mu_{i} + eta_{t} + varepsilon_{it}$$

(7)

Equation (5) is the spatial autoregressive *β* convergence model (SAR), and Eqs. (6) and (7) are the spatial error *β* convergence model (SEM) and the spatial Dubin *β* convergence model (SDM), respectively. In these equations, *ρ* represents the spatial effect coefficient of the explained variable, reflecting the influence of the explained variable in the neighboring cities; *λ* represents the spatial effect coefficient of the error term, reflecting the random shock; *θ* represents the spatial effect coefficient of the explanatory variable, reflecting the influence of the explanatory variable of neighboring cities; and *W*_{ij} represents the spatial weight matrix. In this paper, we use the inverse distance weight matrix.

After estimating the three above spatial panel models, following Alhorst [59] and Han [60], we test and select appropriate spatial econometric models. First, the Lagrange multiplier (LM) test is performed based on the traditional panel model. If the LM test is significant, the SDM model will be used. Then, the likelihood ratio (LR) test and Hausman test are used to determine individual/time effects and fixed/random effects. Subsequently, the Wald and LR tests are used to determine whether the SDM model can be degraded to the SEM or SAR model. If the degraded result is inconsistent with the LM test result, the SDM model will be selected. In the conditional *β* convergence estimation, we also use the same strategy to determine the specific form of the spatial econometric model.

### Data and variables

According to the State Council (2019) [26], the Yangtze River Delta (YRD) region includes Shanghai (SH), Jiangsu (JS), Zhejiang (ZJ), and Anhui (AH), consisting of 41 cities at and above the prefectural level. The main variables and data used in this study are reported below.

Because healthcare resources include human, material, and financial resources [11], two indicators are selected for each of the above three aspects according to the availability of data, and the entropy weight method is used to construct the HRS index [22, 61]. The data source and weight are reported in Table 1. The HRS results of the 41 YRD cities from 2007 to 2019 are reported in Table 11 in Appendix. It can be seen that the HRS of the YRD and sub-provincial regions has shown a continuous upward trend with certain differences among the provinces.

In order to estimate conditional *β* convergence, it is necessary to control the heterogeneity factors across cities. According to the availability of data, we include the following control variables in the model: city size (SIZE), economic development level (GDP), urbanization level (UR), government financial capacity (GFC), and government financial self-sufficiency rate (GFS) [62,63,64]. The definitions and data sources of these variables are reported in Table 2. Logarithm processing is used to eliminate the dimensional difference.

Because the existence of medical universities/colleges may affect the HRS of that specific city, we differentiate the city subsamples according to the existence of medical universities/schools and medical graduate programs in the starting year of the sample. The medical education data of different cities are reported in Table 12 in Appendix.

We also incorporated the main healthcare reforms into the empirical model. The policy variables include the descending healthcare resources reform (DHR) and comprehensive medical reform (CMR). The descending healthcare resources reform is to encourage the balanced distribution by building medical treatment combinations and increasing government investment [65]. Different cities usually experienced different timing of reforms. We take the time of first document issued by the local health departments (authorities) or the time of the first medical treatment combination being organized as the reform time point (see Table 12 in Appendix). We then define one policy dummy DHR, setting the after-reform period as 1 and the pre-reform period as 0.

The comprehensive medical reform tries to reduce diagnosis and treatment costs and realize the goal of hierarchical diagnosis and treatment. In February 2015, Jiangsu and Anhui were identified as pilot provinces for comprehensive medical reform, followed by Zhejiang and Shanghai in May 2016. The date of publication by local health authorities or the official launch date of the reform is used as the time point of reform, so we have the dummy CMR. It should be noted that, if the above reform is launched before July 1 of the current year, the value of 1 is assigned for the current year, and if the reform is launched after July 1 of the current year, the value of 1 will be assigned from the next year.