The proposed Hierarchical 2SFCA method in this study is developed based on several important advancements of 2SFCA, including the Modified 2SFCA [31], the Multiple Catchment Sizes 2SFCA [30], and the Gaussian 2SFCA [24]. The formulas of these models are uniformly expressed using the Generalized 2SFCA framework [1].

### Generalized 2SFCA

The original 2SFCA adopts a dichotomous form of distance decay function, which has been considered a major limitation. Several improved forms of 2SFCA have been developed to address this limitation. The Generalized 2SFCA proposed by Wang [1] can uniformly formulate the various forms as follows:

$$ {A}_i={\sum}_j\frac{S_jf\left({d}_{ij}\right)}{\sum_k{P}_kf\left({d}_{kj}\right)} $$

(1)

where *A*_{i} is the accessibility at location *i*, *S*_{j} is the size of supply of facility *j*, *P*_{k} is the demand size at spatial unit *k*, *d*_{ij} (*d*_{kj}) is the travel cost measured by distance, travel time or total cost, and *f* is distance decay function, which can take discrete forms or continuous forms. In this study, the Gaussian function is adopted to model the distance-decay effects as suggested by existing studies [4, 24]. The Gaussian function can be formulated as:

$$ f\left({d}_{ij}\right)=\left\{\begin{array}{c}\frac{e^{-1/2\times {\left({d}_{ij}/{d}_0\right)}^2-}{e}^{-1/2}}{1-{e}^{-1/2}},{d}_{ij}\le {d}_0\\ {}\kern3em 0,\kern1.5em {d}_{ij}>{d}_0\end{array}\right. $$

(2)

where *d*_{0} is the size of catchment area, i.e., the threshold distance or travel time. Only one parameter is needed for the Gaussian function, making it superior to other alternatives due to less subjectivity and uncertainty caused by parameter setting.

### Relative distance Bias and M2SFCA

The 2SFCA family of metrics, including original 2SFCA and most variants (especially those with distance decay function improvements), have a common nature: the weighted sum of accessibility scores at all locations is equal to the total supply [54]. When calculating the accessibility, the resources at each facility are allocated to the demand nodes within its catchment area based on certain weights, which are determined by the population and the relative distances from the residence locations to the facility. Delamater [31] stated that these metrics generate container-like outputs. The units of spatial accessibility estimated by these methods are the average supply (e.g., physicians) per person. This makes the results of 2SFCA metrics easy to interpret and convenient for policy making. However, this nature of 2SFCA also causes a crucial limitation, which can be explicitly demonstrated by using a simple numerical example shown in Fig. 4.

Suppose there are two different systems A and B. In each system, there is one facility (a) and three demand nodes (x, y, z). The supply of the facility and the demand sizes of three demand nodes are all set as 1. Demand nodes x, y, and z are closer to the facility in system A than in system B. In system A, the values of distance decay function (rather than the distance/travel time itself) for x, y, and z are 0.1, 0.2, and 0.3 respectively. In system B, these values are 0.2, 0.4, and 0.6 respectively, two times higher than in system A.

Based on the above settings, the spatial accessibility of each demand node in system A can be calculated as follows:

$$ {\displaystyle \begin{array}{c}\mathrm{accessibility}\ \mathrm{of}\ \mathrm{x}=1\ast 0.1/\left(1\ast 0.1+1\ast 0.2+1\ast 0.3\right)=1/6;\\ {}\mathrm{accessibility}\ \mathrm{of}\ \mathrm{y}=1\ast 0.2/\left(1\ast 0.1+1\ast 0.2+1\ast 0.3\right)=1/3;\\ {}\mathrm{accessibility}\ \mathrm{of}\ \mathrm{z}=1\ast 0.3/\left(1\ast 0.1+1\ast 0.2+1\ast 0.3\right)=1/2.\end{array}} $$

Similarly, the spatial accessibility in system B can also be calculated as follows:

$$ {\displaystyle \begin{array}{c}\mathrm{accessibility}\ \mathrm{of}\ \mathrm{x}=1\ast 0.2/\left(1\ast 0.2+1\ast 0.4+1\ast 0.6\right)=1/6;\\ {}\mathrm{accessibility}\ \mathrm{of}\ \mathrm{y}=1\ast 0.4/\left(1\ast 0.2+1\ast 0.4+1\ast 0.6\right)=1/3;\\ {}\mathrm{accessibility}\ \mathrm{of}\ \mathrm{z}=1\ast 0.6/\left(1\ast 0.2+1\ast 0.4+1\ast 0.6\right)=1/2.\end{array}} $$

As a result, the accessibility scores of three demand nodes estimated by traditional 2SFCA are the same in both systems. According to our common sense as well as the definition of spatial accessibility, the spatial accessibility in system B should be lower than that in system A, because obviously, spatial barriers are stronger in system B, while other components remain the same in both systems. However, the traditional 2SFCA methods family cannot reflect this difference. The reason is that the term of travel distance works as “relative distance” in these methods [31]. What really matters in the assessment of spatial accessibility from a demand node to a facility is the relative value of the distance from the node to the facility compared to the distances for other nodes within the catchment area of the facility, rather than the absolute value of the distance. The resulting bias can be termed as the “relative distance bias”.

Delamater [31] proposed a Modified 2SFCA (M2SFCA) method which incorporates an additional term of absolute distance into the 2SFCA method. By doing so, spatial accessibility will depend on both relative distance and absolute distance between demand nodes and facilities. According to M2SFCA, the spatial accessibility in system B is lower than that in system A as expected.

M2SFCA can be expressed as:

$$ {A}_i=\sum j\frac{S_jf\left({d}_{ij}\right)f\left({d}_{ij}\right)}{\sum_k{P}_kf\left({d}_{kj}\right)} $$

(3)

where the meanings of all variables are the same as Eq. (1). In Eq. (3), the absolute distance is represented by *f*(*d*_{ij}). Given that *f* is a monotone decreasing function, a larger distance *d*_{ij} indicates a smaller value of *f*, and thus lower accessibility. The term \( \frac{f\left({d}_{ij}\right)}{\sum_k{P}_kf\left({d}_{kj}\right)} \) describes the impact of relative distance on accessibility. The numerator denotes the absolute distance effect for demand node *i*, and the denominator denotes the sum of absolute distance effects for all nodes within the catchment area of facility *j*, which are weighted by the amount of demand. Therefore, the whole term indicates the relative value of a node’s distance compared to the distances of other nodes.

### Relative distance Bias in hierarchical systems

Delamater [31] explained that it is assumed the facility system is optimally configured in traditional 2SFCA methods, whereas in reality, the system configuration is sub-optimal in most situations. The only situation where the configuration of facility system is optimal is when all demand nodes are overlapped with facilities, i.e., the distances from demand nodes to facilities are zero. In a sub-optimally configured system, resources of facilities cannot be fully exploited. The larger the separation between demand nodes and facilities, the lower the exploitation rate of resources.

However, Delamater [31] only considered the situation of single-level facilities and attributed the relative distance bias to sub-optimal configuration of facility systems. It will be shown in this study that the relative distance bias also exists in a hierarchical facility system, albeit in a different way. This can also be demonstrated using a numerical example shown in Fig. 5.

This example deals with a bi-level system. In system C, there is one level-2 (the higher level) facility (a) whose supply size is 8. In system D, there are two level-1 facilities (b1, b2) whose supply sizes are both 4. Both systems have three demand nodes (x, y, z) with their demand size as 1. In system C, the values of distance decay function are 0.2, 0.2, and 0.4 for demand nodes x, y, and z, respectively. In system D, these values are 0.1, 0.1, and 0.5 for x, y, and z for facility b1, and 0.2, 0.2, and 0.1 for x, y, and z for facility b2, respectively. According to traditional 2SFCA methods, accessibility in both systems would be the same, with the accessibility of x, y, and z being 2, 2, and 4 respectively. However, in system C, the spatial barriers from x, y, and z to facility a are 0.2, 0.2, and 0.4 respectively. In system D, the average spatial barriers from three nodes to the two facilities are 0.15, 0.15, and 0.3 respectively. That is to say, in system D demanders are closer to facilities on average than in system C. The spatial accessibility in system D should be higher than that in system C with absolute distance effect taken into account.

The above example explicitly reveals the existence of relative distance bias in hierarchical systems. This is embedded in the allocation of resources among different levels of facilities. In a hierarchical system, higher-level facilities often have larger sizes but a smaller number of facilities, and therefore are on average more distant from demanders that are sparsely distributed across the space. As a result, given a certain amount of resources, i.e., physicians in healthcare facilities, if these resources are configured at a lower level, they are supposed to be closer to demanders on average and therefore generate higher accessibility. In a hierarchical system, both the relative distance and absolute distance effects should be considered in measuring spatial accessibility to eliminate the relative distance bias generated by traditional 2SFCA methods. M2SFCA can be utilized to achieve this goal and to generate better estimation of the spatial accessibility to hierarchical facilities.

### Multiple catchment sizes 2SFCA

Another important feature of hierarchical facilities is that facilities at various levels often have different service scopes. Generally, higher-level facilities serve demanders within larger service areas. The variable catchment areas should be incorporated into the measurements of spatial accessibility to hierarchical facilities. Tao et al. [30] developed a Multiple Catchment Sizes 2SFCA method (MC2SFCA), which determines catchment area sizes according to facilities’ sizes. They supposed that larger facilities generally have larger catchment areas that provide services for more demanders. The method was validated by using a case of residential care facilities in Beijing.

This study proposes to modify the MC2SFCA method for measuring the spatial accessibility to hierarchical facilities. When applied to a hierarchical system, the catchment area sizes are determined based on the hierarchies of facilities, rather than their sizes as in the original form. Higher-level facilities should be assigned with larger catchment areas. Furthermore, though no distance decay parameter is specified in Gaussian-based 2SFCA models, the catchment size parameter (*d*_{0} in Eq. (2)) also determines the slope of distance decay function. A larger catchment size means a weaker distance decay effect. In this sense, the distance decay function is flatter for higher-level facilities. These assumptions are in accord with the central place theory and has been verified by the empirical analysis of Jia et al. [50].

### Hierarchical 2SFCA

The proposed Hierarchical 2SFCA (H2SFCA) method is developed by combining Gaussian-based 2SFCA, M2SFCA, and MC2SFCA, which have been described above. These methods are redefined in hierarchical systems, to reflect several important features of hierarchical facilities. The proposed H2SFCA method can be expressed as:

$$ {HA}_i={\sum}_l{\sum}_{j\in \left\{{d}_{ij}\le {D}^l\right\}}\frac{S_j^lf\left({d}_{ij}\right)f\left({d}_{ij}\right)}{\sum_{k\in \left\{{d}_{kj}\le {D}^l\right\}}{P}_kf\left({d}_{ij}\right)} $$

(4)

where *HA*_{i} is the general spatial accessibility at location *i* in a hierarchical system, *l* is the level of facilities in the system, \( {S}_j^l \) is the supply of facility *j* at level *l*, *D*^{l} is the travel time threshold (i.e. catchment area size) for facilities at level *l*. The distance decay function *f* takes a Gaussian function form as follows:

$$ f\left({d}_{ij}\right)=\left\{\begin{array}{c}\frac{e^{-1/2\times {\left({d}_{ij}/{D}^l\right)}^2-}{e}^{-1/2}}{1-{e}^{-1/2}},{d}_{ij}\le {D}^l\\ {}\kern4em 0,\kern1.25em {d}_{ij}>{D}^l\end{array}\right. $$

(5)

where all variables are the same as Eq. (4).

### Estimation of travel time and catchment sizes

This study uses the Driving Searching API provided by Baidu Map to estimate travel time from demand locations to hierarchical facilities. API allows researchers to obtain reliable estimations of travel time by using the data provided by map service providers [55]. Descriptions of Baidu Map API can be found in previous studies [37, 56]. The estimation of travel time is between 10 a.m. and 5 p.m. on work days to avoid the impacts of traffic congestion during rush hour, which is more related with commuting travel rather than travel to hospitals.

Different transport modes should be considered when measuring the travel costs to facilities at different levels in a hierarchical system. In general, most patients travel relatively long distances to reach tertiary and secondary hospitals by using motorized transport. Therefore, the travel times from demand locations to tertiary and secondary hospitals were estimated using the above API approach. Primary facilities (mainly CHSCs) are roughly configured at the community level and are relatively close to patients’ residential locations. The mean distance from each community to the closest primary facility is 0.55 km. 96.4% communities (743 out of 771) are located less than 2 km from the closest primary facility. Therefore, non-motorized transport such as walking and cycling is preferred for visiting primary facilities. In recent years, due to the “boom of shared bikes” such as Mobike and OFO, cycling has become popular again for daily trips in Chinese cities. We consider cycling as the main transport mode for trips to primary facilities. According to a recent survey report [57], the average speed of shared bike riders in Chinese cities is 9 km/h. The distances from demand nodes to primary facilities are estimated based on actual road networks by using ArcGIS 10.5.

The proposed H2SFCA also assumes that facilities at various levels have different service areas (i.e., catchment areas) in a hierarchical system. Referencing existing studies [37, 56], the catchment area size was set as the maximum travel time by public transport from each demand node to the closest facility, so that all demanders at each location can access at least one healthcare facility within the catchment area. By doing so, the catchment area sizes *D*^{l} are set as 25 min, 40 min, and 70 min for facilities at three levels. Based on the above data and settings, the spatial accessibility at each demand location (i.e., communities) to each level of healthcare facilities, as well as the general accessibility, can be estimated. Since the outputs were accessibility scores at discrete locations, the spatial interpolation analysis is further conducted to generate continuous surface of spatial accessibility distribution. Specifically, the Inverse Distance Weighted method provided in ArcGIS 10.5 was applied to do interpolation.