Revealed by the systematic review, data envelopment analysis (DEA) is the most common approach for assessing efficiency in the health care sector (12). It is a non-parametric technique for measuring the relative efficiency of various entities (referred to as decision making units, DMUs) based on multiple inputs and outputs. The DMUs are sets of homogenous samples, which could be the health system (13, 14), hospital (15, 16), or one specific department/service (17, 18). For instance, Gavurova et al. (14) assessed and compared the national health systems’ performance across OECD countries. While, Akkan et al. (18) took the emergency departments in hospital as the DMUs for efficiency assessment by means of DEA.

There are numerous options available for the specific DEA model. The conventional DEA models indicate to Charnes-Cooper-Rhodes (CCR) and Banker-Charnes-Cooper (BCC) models, which remain the predominating choice for health care efficiency (12). Several pitfalls in the conventional models have gradually been noticed, meanwhile, some derived DEA models emerge and are increasingly adopted in researches. Ferreira et al. (19) applied the non-radial directional distance function model to evaluate the performance of public hospitals in Portugal. This approach takes into account both radial and slack improvements, rather than solely focusing on radial improvement as seen in conventional models. In another study, slacks-based-measure (SBM) DEA model was employed to evaluate the efficiency of Spain hospitals (20). Apart from considering the radial and slack improvements, the SBM DEA is capable of dealing with undesirable outputs. However, in most cases, plural DMUs are estimated as technically efficiency, abovementioned models fail to further rank these efficient DMUs. To overcome this limitation, some scholars utilized super-efficiency DEA model to evaluate efficiency of medical health services (21).

Besides, more complexed DEA models have been developed by integrating diverse derivative models, thereby augmenting their advantages. Thereinto, the super SBM DEA model, which combines the merits of both the SBM and super DEA models, has recently become the preferred choice for efficiency evaluation in many fields, including health care (22), environment (23), urban development (24), and logistics industry (25). As for the Chinese health system, Sun et al. conducted a super SBM DEA model to explore and compare the efficiency of health care services across 31 provinces. Similarly, Zhao et al. employed the same approach to measure the performance of primary health care institutions in China. After analyzing the distinctive features and applications of various common DEA models, we determined to construct a super SBM DEA model in our study for measuring the efficiency of specialized public health facilities.

In practice, most researchers concentrate on the efficiency assessment for medical services (12, 26), while attention are comparatively inadequate for preventive public health services. From a global perspective, there are significant geographical disparities in the performance of public health systems, for instance, public health systems in the European region are more efficient than those in the African, South-East Asia, and Western Pacific regions (27). By means of DEA, Mukherjee K et al. found that public health sectors operated with nearly 28% inefficiency in the US (28), and emphasized the significance of internal funding for solving this inefficiency. Similarly, there was 24% inefficiency in the public health sectors of Senegal (29). Besides, a recent study evaluated the performance of the unintentional childhood injury intervention, one public health program in Japan, which presented with increasing trends and shrinking regional disparities (30).

Similarly, the performance of specialized public health facilities, the principal supplier for public health services in China, has been less extensively studied compared to that of hospitals and primary health care institutions (31–37). By employing the traditional DEA model, Dong LM and Zu JT discovered that specialized public facilities demonstrated a lower level of efficiency in contrast with primary health care facilities and hospitals in China from 2009 to 2017 (31). Using the same approach, scholars have reported the inefficiency of 56% and 55.6% in specialized public health facilities in 2015 and 2019 (38, 39), indicative of remarkable underutilization of available public health resources. Besides, some assessed the efficiency of certain specialized public health institution, such as centers for disease control and prevention (40) and maternal-child health institutions (41). So far, only a few studies measured the efficiency of specialized public health facilities in China, and mostly employed the traditional DEA model.

Understanding the determinants of efficiency could provide empirical evidence for policy formulation. Given limited researches exploring the determinants affecting public health efficiency, this subsection provides a review of the methods and findings related to health care system.

It has been extensively discussed how demographic, social, and economic factors affect the efficiency of health systems. For example, urbanization rate was found to be positively associated with technical efficiency of primary health care system in China (42). A higher proportion of vulnerable population, like children or older adults, tended to reduce the efficiency of health care system (2). The relationships between different sources of health spending and health system efficiency have empirically been recognized (43). Besides, the environmental elements are associated with population wellness, which may affect the performance of public health system (44).

Econometric models, including Tobit regression (45, 46), bootstrap truncated regression (47), and multiple linear regression (48), are usually employed to investigate the potential determinants for health system efficiency. However, these models estimate the regression parameters as a whole, failing to distinguish the unique impacts on different sub-regions within a single regression analysis. In a prior Chinese study, it constructed four Tobit regression models to analyze the determinants for rural public health efficiency from perspectives of national, eastern, central and western regions, and remarkable differences were noted in the regression results. Specifically, the economic level was remarkably associated with the rural public health services in eastern and western regions, not the central region (49). This finding confirms our concern that a specific determinant may have distinct impacts in different geographic regions.

The geographically weighted regression (GWR) model enables the parameter estimation from the global and local spatial scales, identifying the significant determinants for each geographic unit. Several scholars constructed GWR model to identify the spatially varied associations between socio-economic factors and health resource allocation in China (50, 51). Therefore, we employed the GWR model in our research to in-depth understand the spatial heterogeneity of factors affecting public health facility efficiency, which also extends the application of GWR approach in health system efficiency.

For efficiency assessment in health system, labor-, capacity- and expense-related indicators are generally selected as inputs, and activity-related indicators constitute the outputs (26). Health technicians which involve doctors, nurses and pharmacists, make up nearly 80% of total personnel and handle the majority of services in specialized public health facilities (38, 39). For this reason, the number of health technicians per thousand people is chosen as the proxy for labor input. Besides, we identified the number of beds per thousand people in professional public health services as the capacity-related input (52, 53). Also, we noted that the expense-related variables had been used in efficiency evaluation for hospitals, primary health institutions and other health services (12, 26, 54, 55), therefore, we involved the operating cost per capita of specialized public health facilities as another input indicator in our study.

In the context of China, specialized public health facilities are primarily responsible for the prevention and management of infectious and chronic diseases, maternity and child hygiene, epidemic surveillance and health education (8). In correspondence to their primary mission, several outputs have been frequently chosen for measuring the efficiency of public health services, including infectious incidence rates of category A and B, maternal health management rate, under-three child health management rate and the number of health education (38, 39). However, the number of health education is more like an intermediate outcome for the former three indicators, which we thought is unreasonable to be listed as one output. Therefore, we confirmed three output variables for specialized public health facilities, namely infectious incidence rate of category A and B, maternal health management rate and under-three child health management rate, which comply with the officially regulated major goals of such agency (7, 56).

Altogether, we included three input and three output variables. And 31 provinces were treated as the DMUs in our study (Figure 2). Input and output variables were retrieved from National Health Financial Annual Report (each specialized public health facilities is required to annually submit key statistics to an electronic system, and we retrieved operating cost per capita of specialized public health facilities from this system) and China Health Statistics Yearbook (57). The detailed statistics of indicators are depicted in Table 1.

Since few studies investigate the influential factors for the performance of public health services, we identified several independent variables associated with health system efficiency (2, 42–44). The independent variables were collected from China Statistical Yearbook (58) and China Statistical Yearbook on Environment (59). Summary statistic for all independent variables is presented in Table 1.

The DEA model has been widely used in various fields since its initial proposal in 1978 (2). The underlying principle of DEA is to construct an optimal production frontier through a linear programming approach. Unlike parametric methods, DEA does not require the pre-specification of any functional form. Besides, it can tackle with multiple inputs and outputs, which also makes DEA stand out among a variety of efficiency evaluation methods. As mentioned before, the conventional CCR and BCC models only assess the radial efficiency of DMUs, which means that the inefficient DMUs is assumed to improve their performance through proportional reduction in the inputs or proportional increase in the outputs. It disregards the slack improvements existing in input or output variables, which possibly leads to an overestimation of the actual efficiency. Moreover, it is quite difficult to control the radial improvements in input and output variables in most practices. Therefore, the non-radial DEA models has been increasingly applied in researches of efficiency evaluation, which could directly handle the potential slack improvement and avoid the radial deviations. Apart from this, the SBM DEA model, belonging to the non-radial model, has been empirically proved to exhibit excellent performance in dealing with undesirable outputs and applied into a mass of academic researches (22–25, 60).

According to the principle of conventional DEA model, when one DMU states in the optimal production frontier, indicating achieving fully efficient with the maximum value of 1. In most cases, plural DMUs have the “efficient status.” Thus, it is meaningful to further distinguish the efficiency among these efficient DMUs. In 2002, Tone proposed an SBM of super-efficiency (hereafter referred to super SBM DEA model), which is capable of ranking the efficient DMUs (61). The infectious incidence rates of category A and B is an undesirable output variable, therefore the super SBM DEA model is selected in our research to handle the issues of undesirable output, slack improvements and ranking efficient DMUs.

Let us assume n DMUs in our research. Each DMU produces s1 desirable outputs and s2 undesirable outputs by means of m input variables. According to Tone (60), the SBM model with can be specified as follows in this case:

Subject to

Where λ denotes the weight vector; S−, Sg, and Sb indicated the slacks of inputs, desirable outputs and undesirable outputs, respectively. And ρk is the technical efficiency of public health service of province k ranging from 0 to 1. If ρk equals to 1, the kth DMU achieves the efficient status. Then, the super-SBM DEA model is applied to these efficient DMUs to further rank their efficiency values. The formula is specified as follows (61):

Subject to

Where ρk∗ indicates the super efficiency of the kth DMU, whose value equals or exceed 1. And the others have the same meaning with abovementioned SBM DEA model with undesirable outputs.

Technical efficiency in DEA analysis refers to “the capacity or even the ability to make the most appropriate use of the resources available to achieve an intended result” (62). Specifically, the technical efficiency is measured under the assumption of constant return to scale. However, not every sample would be stated in the aforementioned status, and the technical efficiency would be influenced by sample’s scale. Hence, the technical efficiency is decomposed into pure technical and scale efficiency. Pure technical efficiency is correlated with managerial factors or technology level (e.g., advanced technology and equipment in health institutions), while scale efficiency implies whether the operation size of sample is optimal, too large or too small, which corresponds to constant-, decreasing- or increasing returns to scale, respectively (63). The relations between three types of efficiency can be specified as follows:

And all the efficiency analysis was conducted in Max DEA software.

Global Moran’s I index was applied to explore the spatial patterns of estimated technical efficiency of specialized public health facilities. The formula of Global Moran’s I was specified as follows (64):

Where, n denotes research samples (31 provinces), Wij represents the matrix of spatial weights and inverse distance weight matrix is adopted in this study to comply with Tobler’s First Law of Geography, Xi and Xj indicates the technical efficiency of province i and province j, respectively.

The value of I ranges from −1 to 1. If the value is equal to 0, no spatial autocorrelation exists; otherwise, positive (I > 0) or negative (I < 0) spatial autocorrelation exists. This procedure was performed in ArcGIS version 10.7.

To choose significant influential factors, multi-variate linear stepwise regression was constructed based on the Z-scores of the dependent variable (technical efficiency obtained from super-SBM DEA model), and all 17 independent variables in Table 1. These analyses were carried out using SPSS version 25.

Although the traditional econometric regression model is usually used to capture the associations between dependent and independent variables. It takes no account of the spatial features of samples. In most cases, each observation contains some geospatial information; for instance, the 31 provinces in our research inherently carry geographical information, which is reflected in the varying distance between the provinces. The ignorance of spatial effects might lead to bias and inconsistency for the parameter estimation. Under this condition, several spatial econometric models, including spatial error model, spatial lag model and spatial Durbin model, have attracted intensive attention with the advantages of handling the spatial interdependency in error terms, dependent or independent variables. However, these spatial econometric models also belong to a global analysis, in another word, a unified parameter is developed for each explanatory variable (65). And the problem of spatial heterogeneity, referring to uneven distribution of measurements over space, still remained unsolved.

The GWR, proposed by Fotheringham, is a statistical technique to discover the spatially non-stationary relationships between dependent and independent variables (66). It embeds geographic characteristics into regression parameters, allowing for the estimation of a unique set of regression coefficients for every geographical unit. This distinctive capability enables GWR to capture local variations in influential factors, thereby assisting local governments in implementing more targeted and effective interventions, as opposed to one-size-fits-all approaches. Consequently, the GWR model was employed in our research to comprehensively identify the spatial heterogeneity between the explanatory variables and the estimated efficiency.

In practice, when spatial autocorrelation is detected, which means the occurrence of spatial non-stationarity, the GWR should be considered (67). The basic structure of GWR is (64):

where, yi is the technical efficiency of province i, xik represents the k^th independent variable for province i, μivi describes the spatial coordinates of province i, β0μivi and βkμivi denote the regression coefficients, and εi is the error term.

The spatially varied parameters are estimated relying on the spatial weight matrix between data points. In this study, corrected Akaike information criterion (AICc) is used to determine the optimal bandwidth and Gaussian kernel function is constructed to calculate the spatial weight matrix. The GWR analysis and mapping of coefficients were performed in ArcGIS version 10.7.

As estimated, the average technical efficiency, pure technical efficiency and scale efficiency scores of specialized public health facilities were 0.6569, 0.7336 and 0.9206, respectively. Regarding the temporal trend, technical efficiency kept decreasing from 0.6889 in 2017 to 0.6238 in 2019, which was primarily attributed to a substantial reduction in the scale efficiency (Table 2).

There exist obvious geographic variations in the efficiency of specialized public health services. The technical efficiency and pure technical efficiency were the highest in the eastern region, followed by the central and western regions in all research years. Remarkably, the pure technical efficiency of the eastern region was nearly 1.5 to 2 times higher than that of the central and western regions. Since the elite public health resources, such as highly skilled professionals, advanced information systems, and cutting-edge management and technologies, are predominantly clustered in the eastern regions. Contrarily, the western region possessed a relatively higher scale efficiency, and experienced less reduction in the scale efficiency over years (Table 2).

The findings at the provincial level provided more insights into geographic variations (Figure 3; Supplementary Table S1). Merely one-quarter of the provinces surpassed the average score in terms of technical and pure technical efficiency, with the majority of them situated along the eastern coast, such as Jiangsu, Shanghai and Tianjin. This precisely unveiled the pressing issues in management and technology within Chinese public health system. Conversely, as regards the scale efficiency, a small fraction of provinces fell below the average score, particularly some economically developed provinces like Shanghai, Zhejiang and Beijing. This indicates that, for prosperous provinces, the emphasis of local governments should shift from continually expanding investment to adjusting the structure of existing resources, such as optimizing personnel quality and expenditure structure.

Based on the stepwise regression, four variables were selected with a 95% confidence interval, which were ratio of older population, sex ratio, proportion of public health expenditure and total volume of Nitrogen Oxides (NO_x) emission. Generally, a variance inflation factor (VIF) larger than 10 implies multi-collinearity. Since the largest VIF was merely 1.283, no multi-collinearity existed in the selected independent variables. The details for stepwise regression are depicted in Supplementary Table S2.

The GWR model is appropriate only when spatial correlation exists in the dependent variables. The global Moran’s I index revealed significant spatial correlations in the technical efficiency (Supplementary Table S3), suggesting the suitability and necessity of GWR.

As a result, the GWR regression was constructed between the technical efficiency and four independent variables. The details of the GWR regression are described in Supplementary Table S4 and Supplementary Figure S1. After the GWR regression, Global Moran’s I of the standard residual became insignificant, reaffirming the reliability of GWR regression (Supplementary Table S5).

According to the GWR regression, the proportion of public health expenditure and ratio of older population had a relatively stronger impact on public health services efficiency as a whole (Supplementary Table S4). During the research period, the proportion of public health expenditure negatively impacted public health efficiency. Its impact slightly decreased over time, whose absolute median value ranged from 0.4524 to 0.5480. A positive association was observed between ratio of older population and public health efficiency, with median coefficients changing from 0.4257 to 0.4793. Similarly, sex ratio was found to positively impact public health efficiency, whose impact reached the largest with a median coefficient of 0.3441 in 2018. The total volume of NO_x emission was the least influential factor with the median coefficient value from 0.1691 to 0.2773.

Subsequently, we displayed the coefficients into maps to explicitly exhibit the impacts of independent variables at provincial level.

For the ratio of older population, its coefficient showed a descending trend from northeast to southwest regions. Specifically, ratio of older population had the highest impact on efficiency in Heilongjiang, Jilin and Liaoning all the time, followed by some other eastern provinces such as Beijing, Tianjin, Hebei and Inner Mongolia. Whereas the weakest association was found in provinces located in the west and the south regions including Xinjiang, Tibet and Hainan (Figure 4).

The positive associations between sex ratio and efficiency gradually weakened from the north to the south. The coefficient of sex ratio was relatively higher in Inner Mongolia, Xinjiang, Gansu and Heilongjiang with a value of around 0.4; on the other side, the lowest coefficient appeared in some southern provinces especially Hainan, Guangdong and Guangxi with an approximate value of 0.2 (Figure 5).

Different from the other three variables, a stronger negative association was identified between the proportion of public health expenditure and public health efficiency. No obvious geographical regularity was found in its coefficient distribution. In 2017 and 2018, the highest absolute coefficient primarily concentrated in central regions such as Shaanxi, Hubei, Henan, and Hunan, while it changed to eastern coastal provinces such as Fujian, Zhejiang, Shanghai, and Jiangxi in 2019. Besides, a comparatively weaker associations were found in provinces including Xinjiang, Tibet, Qinghai and Gansu, whose absolute coefficients ranged around 0.25–0.4 all years (Figure 6).

The association between total volume of NO_x emission and technical efficiency was the least among four independents. Notably, its impact gradually dampened from southeast to northwest. Higher coefficient was generally observed in provinces such as Guangdong, Hainan, Fujian and Jiangxi with a value of 0.2–0.3, whereas provinces with lower coefficients mostly distributed in the northwest region including Xinjiang, Gansu, Tibet and Qinghai (Figure 7).