Scale-Dependent Spatial Clustering of Hospital Operational Performance and Patient Safety Outcomes in East Java, Indonesia: A Spatial Autocorrelation Analysis

1. Introduction

Hospital performance and patient safety are central concerns for health systems worldwide. According to the Global Patient Safety Report 2024 by the World Health Organization (WHO), unsafe care causes more than 3 million deaths annually worldwide [1]. In low-and middle-income countries alone, 134 million adverse events occur annually in hospitals, leading to an estimated 2.6 million deaths. While considerable attention has been paid to improving patient safety in high-income settings, critical gaps exist in hospital performance and patient safety culture in low-and middle-income countries due to limited capacity and health budget constraints [2].

Indonesia is the world’s fourth most populous country, and has a decentralized healthcare system [3]. According to Indonesia Health Profile 2024, there are 3288 hospital facilities nationwide, each with distinct ownership, administrative class, hospital type, and accreditation status,- all of which may influence care quality in distinct ways [4]. Approximately 95.3% of hospitals are accredited [4], though average compliance with national quality indicator reporting remains around 69% [5]. Hence, accreditation status alone does not guarantee full compliance with operational performance or patient safety outcomes. Analyses of reporting compliance have further revealed significant variation across provinces and hospital classes [6]. Additionally, there is evidence of variation in healthcare outcomes at the district level [7,8]. Such heterogeneity raises important questions about whether these disparities in hospital performance and patient safety follow a geographic pattern at the hospital or district scale.

Geospatial methods offer a powerful approach to answering this question. Spatial autocorrelation analysis, including global and local Moran’s I, enables the detection of non-random geographic clustering in health outcomes, thereby identifying hotspots of elevated risk or concentrated underperformance [9,10]. Without geospatial analysis, such hot and cold spots cannot be specifically visualized and remain hidden within national data [11]. Beyond spatial autocorrelation, geographic information systems (GIS) have been applied extensively to investigate healthcare systems from various perspectives. For instance, Cuadros, et al. [12] used GIS for early detection and management of viral outbreaks. Similarly, Weiland, et al. [13] identified geographical inequalities that led to avoidable infant mortality. Furthermore, Moragues, et al. [14] investigated the accessibility of public hospital services in Mallorca using GIS. Recent applications of GIS for public health also include research assessing spatial accessibility and equity of public hospitals for older adults in Shanghai [15]. Notably, GIS-based approaches were adopted widely and rapidly in response to COVID-19, underscoring the usefulness of spatial analysis for public health [16,17]. Additional evidence supports the application of GIS in public health research across a wide range of countries including Japan [18], Pakistan [19], and Canada [20]. However, despite this growing adoption, its application to hospital performance and patient safety outcomes is lacking in the Indonesian context.

This study addresses this gap by conducting a comprehensive spatial analysis of 428 hospitals across East Java province. Standardized composite indices of operational performance and patient safety outcomes are constructed from routinely collected indicators including bed occupancy rate (BOR), average length of stay (ALOS), bed turnover (BTO), gross death rate (GDR), and net death rate (NDR). These indices are used to examine whether hospital performance and patient safety outcomes are spatially clustered at the hospital or district scale. The findings have direct relevance for health system planners seeking to identify geographical disparities, strengthen quality governance, and prioritize regulatory efforts where risk is greatest.

2. Materials and Methods

2.1. Study Area

This research was conducted in the East Java province, one of Indonesia’s 38 provinces. It is among the most populous provinces and has the highest concentration of hospitals as well as specialist and advanced services [21]. The hospital network in East Java is structurally diverse, spanning urban centers and rural district while covering a broad spectrum of ownership types, accreditation levels, and class categories. These characteristics make East Java an ideal place to investigate whether geographic inequalities manifest in hospital operational performance and patient safety outcomes.

2.2. Data Acquisition and Explaination of Variables

Data for the year 2024 were obtained from the East Java Provincial Health Office. These data are collected as part of monthly reports submitted by hospitals to the provincial authorities. The dataset comprised 445 hospitals, each listed with hospital identification variables, institutional attributes, operational performance indicators, and patient safety indicators. All indicators and their definitions are listed in Table 1.

2.3. Data Cleaning and Preprocessing

The data were stored in a spreadsheet. Using filtering and visualization, nine hospitals were found to have missing data relating for operational performance and patient safety indicators; these hospitals were excluded from further analysis. Additionally, eight hospitals were identified as “newly opened” hospitals using the remarks column. As these hospitals were new, their value against indicators were zero and were excluded from further analysis. Ultimately, 428 hospitals remained and were included in the final analysis.

Further, the 17 types of ownership were grouped into 6 analytically meaningful ownership groups: private sector, local government, central government, faith-based organizations, social/non-governmental organizations (NGOs), and military and police.

Finally, the latitude and longitude of all hospitals were recorded through geocoding process, which involves converting location descriptions - such as place name, street address, and postal code - into geographic coordinates [24]. Using hospital names and district information, each hospital was located via Google Maps and coordinates were recorded in the spreadsheet. This step was crucial for enabling data to be represented in a GIS environment for spatial analysis.

2.4. Transformation and Standardization of Indicators

The ALOS, GDR, and NDR are indicators where higher values represent poorer outcomes. Conversely, higher values for BOR and BTO indicate better outcomes. To ensure consistency and comparability, the values for ALOS, GDR, and NDR were inverted so that higher values indicate better outcomes. Additionally, BOR was transformed based on its deviation from the optimal threshold (75%).

In the next step, all five indicators were standardized using z-score normalization [25]. This procedure was performed separately for each hospital type (general, maternity, and specialist) to account for structural and functional differences across hospital types, differing significantly in terms of patient case mix, service complexity, and operational characteristics. Given that these differences can lead to systematic variation in performance and safety indicators, applying a single global standardization across all hospitals could bias the results by unfairly comparing inherently different hospital types. Therefore, type-specific standardization ensures that each hospital is evaluated relative to its peer group, enabling more meaningful comparison and reducing structural bias in the final composite indices. The transformation is conducted as per Equation 1:

$$Z_{it}={\displaystyle \frac{X_{it}-{\mu }_t}{{\sigma }_t}}$$

(1)

where$X_{it}$= value of indicator for hospital $i$in type $t$,${\mu }_t$= mean of indicator within hospital type $t$, and${\sigma }_t$= standard deviation within hospital type $t$.

Figure 1 offers an overview of the main steps of data processing and analysis.

2.5. Construction of Composites

Rather than relying on individual indicators, composite indices offer an excellent approach to reliably rank hospitals within certain contexts (e.g., operational performance or safety risk) [26]. The five indicators in our dataset correspond to two different themes, i.e., hospitals’ operational performance and patient safety outcomes. Thus, two separate composite indices were developed to summarize individual indicators into these themes while preserving their conceptual distinction. These indices were created following the established methodologies for indicator normalization and aggregation outlined in the OECD’s Handbook on Composite Indicators [27] and widely applied in health system performance assessment frameworks [28].

2.5.1. Composite Operational Index

The composite operational index (COI) was constructed to capture hospital efficiency and service delivery performance using standardized values of BOR, ALOS, and BTO. It was calculated as the average of the standardized values (Equation 2), providing a unified measure of operational performance. Higher values indicate better efficiency and utilization of hospital resources.

$$CO{I}_i={\displaystyle \frac{Z_{BOR,i}+Z_{BTO,i}+Z_{ALOS,i}}{3}}$$

(2)

where $Z_{BOR,i}$, $Z_{BTO,i}$, and $Z_{ALOS,i}$ represent standardized values of BOR, BTO, and ALOS, respectively.

2.5.2. Patient Safety Index

The patient safety index (PSI) was constructed to represent hospital safety outcomes using directionally adjusted and standardized values of GDR and NDR. In the final PSI, higher values indicate better patient safety. The composite index was computed as the average of the standardized indicators (Equation 3), providing an integrated measure of hospital-level safety outcomes.

$$PS{I}_i={\displaystyle \frac{Z_{GDR,i}+Z_{NDR,i}}{2}}$$

(3)

where $Z_{GDR,i}$and $Z_{NDR,i}$ denote standardized values of directionally adjusted GDR and NDR.

2.6. Descriptive Statistics

Descriptive statistics were computed to summarize the distribution of hospital characteristics and the five indicators. Continuous variables were summarized using means and standard deviations, whereas categorical variables were presented as frequencies and proportions. The distribution of hospitals across ownership groups, hospital types, and class categories was examined to provide an overview of the study sample. Additionally, the distribution of standardized indicators and composite indices was assessed to identify potential outliers and verify comparability across hospital types following standardization.

2.7. Spatial Distribution of Hospitals

To visualize the distribution of high and low performing hospitals with respect to operational performance and patient safety, COI and PSI were classified into four categories - high (≥ 1), moderately high (0 - < 1), moderately low (> -1 - < 0), and low (≤ −1) -following a quartile-based classification approach commonly utilized in studies employing composite indices [29,30].

Figure 2 outlines the workflow of spatial autocorrelation analysis. To evaluate the spatial dependence of hospital operational performance and patient safety indicators, global spatial autocorrelation was quantified using Moran’s I statistic. This metric tests whether observed values are spatially clustered, dispersed, or randomly distributed across geographic space [31,32]. If spatial clustering is evident in global Moran’s I, then local Moran’s I (Anselin local Moran’s I) can further reveal the underlying spatial heterogeneity. Specifically, local Moran’s I—also called cluster and outlier analysis—identifies concentrations of high and low values and spatial outliers [33,34].

Initially, point-based global Moran’s I analysis was conducted at the hospital scale separately for COI and PSI values. To further disaggregate the spatial dynamics of mortality, global Moran’s I was also independently calculated for standardized GDR and NDR. If evidence of clustering was found in any variable, local Moran’s I analysis was conducted.

Additionally to hospital scale, global spatial autocorrelation was subsequently evaluated at district level to assess broader regional trends. Using a spatial join procedure, district-level mean composite index values were aggregated from the hospital data. If any evidence of clustering emerged, local Moran’s I analysis was also conducted to district-level composite index averages. Furthermore, to analyze scale-dependent patterns of mortality, local Moran’s I was analyzed individually for standardized GDR and NDR at district scale. The analysis classified districts into four categories: high–high (HH), low–low (LL), high–low (HL), and low–high (LH), representing clusters of similar values and spatial outliers.

3. Results

3.1. Descriptive Statistics and Spatial Distribution of Hospitals

A total of 428 hospitals from across East Java province were included in the analysis. As illustrated in Figure 3(a), general hospitals were by far the predominant hospital type in the sample (n = 346, 80.8%), followed by maternity hospitals (n = 63, 14.7%); specialist hospitals represented only a small proportion (n = 19, 4.4%).

Regarding hospital classification (Figure 3[b]), most facilities were categorized as Class C (n = 224, 52.3%), followed by Class D (n = 134, 31.3%), and Class B (n = 63, 14.7%); Class A hospitals constituted a very limited share (n = 7, 1.6%). According to Indonesian regulations, hospital classification is partly determined by bed capacity; Class A hospitals have at least 250 beds, Class B hospitals have 200-249 beds, Class C hospitals have 100-199 beds, and Class D Hospitals have a minimum of 50 beds.

In terms of ownership structure (Figure 3[c]), most hospitals were run by the private sector (n = 260, 60.7%), with a substantial contribution from local government hospitals (n = 86, 20.1%). Other ownership categories included faith-based organizations (n = 27, 6.3%), military and police institutions (n = 26, 6.1%), and social/NGO-owned hospitals (n = 22, 5.1%), while central government hospitals accounted for only a marginal proportion (n = 7, 1.6%).

As shown in Figure 3(d), the vast majority of hospitals had achieved plenary accreditation or Paripurna status (n = 370, 86.4%), representing the highest level of accreditation. This indicates a high level of compliance with national healthcare standards. In contrast, main accreditation or Utama status—reflecting a moderate level of compliance—was observed in 10.7% (n = 46 hospitals); only 2.8% (n = 12 hospitals) received intermediate accreditation (Madya status).

The spatial distribution of hospitals across districts is shown in Figure 4. Only one district, Kota Surabaya (n = 62), was found in the highest range of hospital concentration (>30–65 hospitals). This is because Surabaya is the capital and largest city of East Java province and the second-largest city in Indonesia, after Jakarta. This was followed by districts in the 20–30 category, including Sidoarjo (n = 30), Kota Malang (n = 28), and Malang (n = 25). A moderate-to-high distribution (>15–20 hospitals) was observed in Gresik (n = 19) and Lamongan (n = 18), while districts such as Jombang (n = 15) and Jember (n = 14) fell within the 10–15 range. The majority of districts were concentrated in the 5–10 category, whereas several districts, including Bondowoso, Nganjuk, and Trenggalek, were classified in the lowest range (3-4).

3.2. Spatial Distribution of Composite Operational and Patient Safety Indexes

The spatial distribution of COI across districts reveals a clear dominance of moderately performing hospitals, with the majority of facilities clustered within the moderately high and moderately low categories (Figure 5). Across most districts, hospitals in the moderately high category consistently outnumber those in the high-performance group. Thus, while operational performance is generally acceptable, only a limited proportion of hospitals achieve optimal efficiency levels. For instance, districts such as Surabaya (n = 21), Sidoarjo (n = 14), and Malang (n = 14) exhibited relatively higher concentrations of moderately high-performing hospitals, reflecting stronger operational capacity in these urbanized and resource-rich areas. However, even within these districts, a substantial number of hospitals remained in the moderately low category, highlighting notable intra-district variability in performance.

Additionally, the high-performance category was sparse across the province, with only a few districts such as Bojonegoro (n = 2), Kota Kediri (n = 2), and Kota Malang (n = 2) exhibiting more than one high-performing hospital; the remaining districts recorded either a single or no hospitals in this category. Conversely, the presence of low-performing hospitals was spatially uneven, with the higher concentration observed in Surabaya (n = 9) and smaller counts in districts such as Bojonegoro, Pamekasan, and Ponorogo. Overall, the spatial pattern indicates that hospital performance across East Java is concentrated at intermediate levels, with limited representation at the extremes, highlighting both the absence of widespread high-performance clusters and the persistence of localized pockets of lower operational efficiency.

Similar to COI, the spatial distribution of the PSI index, where higher values indicate lower mortality rates and better patient safety performance, reveals predominantly intermediate safety levels across East Java districts (Figure 6). Most hospitals were clustered within the moderately high and moderately low categories, indicating that patient safety outcomes are generally acceptable but not consistently optimal. Districts such as Surabaya (moderately high: n = 31) and Sidoarjo (moderately high: n = 13) exhibited relatively strong safety profiles, reflecting better clinical performance in highly urbanized and resource rich settings. Similarly, Malang (high: n = 7; moderately high: n = 7) and Gresik (high: n = 4; moderately high: n = 9) demonstrated comparatively favorable distributions, suggesting concentration of hospitals with stronger outcomes.

Despite these encouraging patterns, hospitals classified within the high-PSI category remained relatively scarce across the province. Notable concentrations were observed only in a limited number of districts, including Malang (n = 7), Sidoarjo (n = 7), and Surabaya (n = 5). In contrast, hospitals with lower-PSI category was more widely dispersed, with the largest number located in Surabaya (n = 8) and Malang (n = 5), as well as in several other districts with smaller clusters. Additionally, a substantial number of hospitals fell within the moderately low PSI category in districts such as Malang (n = 11) and Surabaya (n = 18), indicating substantial variations in patient safety outcomes within the same districts. Overall, the findings suggest that while extremely poor safety outcomes are relatively uncommon, the absence of a strong concentration of high-performing hospitals and the persistence of intra-district disparities indicate significant opportunities for improvement. Strengthening quality and safety initiatives, particularly in hospitals performing within the moderately low category, may help reduce geographic variation and support more consistently high patient safety standards across the province (Figure 6).

3.3. Outcomes of Global Moran’s I

Table 2 summarizes global Moran’s I results at hospital and district level. The spatial analysis showed different patterns between hospital-level and district-level indicators. At the hospital level, neither the COI (Moran’s I = -0.005, p = 0.905) nor the PSI (Moran’s I = -0.025, p = 0.218) demonstrated significant spatial autocorrelation, indicating that hospital performance and mortality outcomes were not geographically clustered or relatively randomly distributed across East Java. The low Moran’s I value indicates that variations in patient safety outcomes are spatially random across individual hospitals. These findings imply that geographic proximity does not play a dominant role in shaping patient safety outcomes at the micro (hospital-level) scale.

In contrast, at the district level, PSI showed significant positive spatial autocorrelation (Moran’s I = 0.348, p = 0.004), indicating that districts with similar patient safety outcomes tended to be located near one another. District with similar poor outcomes were more likely to be surrounded by districts with similar outcomes. This findings suggests that patient safety outcomes, in this case is the mortality pattern, maybe influenced by broader contextual factors operating at district level, such as healthcare resources, referral systems, worforce availability, and socioeconomic conditions.

Because PSI was derived from two mortality indicators, GDR (Moran’s I = 0.353, p = 0.004), and ND (Moran’s I = 0.342, p = 0.005), these indicators were analyszed separately to assess the robustness of the findings. Both GDR and NDR also demonstrated significant spatial clustering. The consistency of these results suggest that the observed spatial pattern in PSI is primarily driven by geopgraphic differences in mortality outcomes across districts.

3.4. Local Spatial Autocorrelation of Individual Patient Safety Indicators

Building upon the global spatial autocorrelation findings, which demonstrated significant district-level clustering for standardized GDR and NDR, Local Moran’s I analysis was conducted to identify the specific location of spatial clusters in district-level GDR and NDR. This analysis provides a more detailed understanding of localized patterns in PSI outcomes across East Java.

For both GDR and NDR, significant low-low (LL) clusters were identified in Mojokerto and Bangkalan (Figure 7 and Figure 8), indicating districts with consistently low mortality rates that were surrounded by neighboring districts exhibiting similarly favorable outcomes. Additionally, a low-high (LH) spatial outlier was observed in Magetan, representing a district with relatively low mortality surrounded by relatively higher-risk neighboring districts. High-high (HH), representing district with elevated mortality rates surrounded by neighboring districts with similarity high mortality, were also detected for both indicators, although their exact locations varied slightly. For GDR, HH clusters were observed in Bondowoso and Jember, whereas for NDR, they were identified in Madiun and Bojonegoro. Importantly, these HH clusters were in adjacent districts, suggesting the presence of a regional hotspot of elevated morality as compared to surrounding regions.

Overall, the appearance of high-risk clusters across both indicators, together with the spatial proximity of low-risk areas, indicates robust and spatially coherent patterns in PSI outcomes at the district level, reflecting underlying regional disparities.

4. Discussion

4.1. Spatial Patterns and Governance Implications

The absence of spatial clustering in hospital operational performance and patient safety indicators at the individual facility level suggests that geographic proximity between hospitals does not explain variations in hospital outcomes. In other words, hospitals with poor performance are not necessarily located near other poorly performing hospitals, and high-performing facilities are not consistently concentrated within specific geographic areas. This pattern indicates that outcomes are driven by factors other than geographical attributes, such as leadership effectiveness, organizational culture, workforce capacity, and institutional governance. In a structurally heterogeneous hospital system; where private, public, military, faith-based, and specialist facilities coexist within the same geographic space - such institutional diversity likely suppresses any detectable spatial signal at the point level [33].

However, the emergence of significant clustering at the district level reveals a different dimension of the health system performance. When hospital-level variation is aggregated, spatial patterns become apparent, suggesting that hospitals within the same districts are influenced by common structural and environmental conditions. These conditions may include governance capacity, resource availability, referral system effectiveness, and broader population health characteristics. This scale-dependent pattern is consistent with spatial epidemiological theory, which recognizes that observed geographic patterns depend on the chosen unit of analysis [35,36]. It also suggests that districts function as governance units that shape a common operating environment for hospitals, ultimately influencing patient safety outcomes.

This interpretation aligns with Indonesia’s decentralized health system in which district governments are responsible for the management of health services. As confirmed by previous studies [8,37,38], the variation in service delivery can be caused by several interrelated factors including health workforce availability, financing, and governance. The challenge of district-level variation in health services has also been documented in Pakistan [39]. Furthermore, the identification of high-risk clusters; where elevated patient safety risk is concentrated across multiple facilities which points to systemic, district-level challenges rather than isolated hospital failures. For policymakers, this distinction is critical. If patient safety outcomes are shaped by shared-district level conditions, interventions focusing exclusively on individual hospitals may have limited effectiveness in addressing the underlying drivers of risk. Instead, geographically targeted governance strategies, prioritizing high-risk districts, may offer a more effective approach to improving patient safety and advancing SDG targets related to health system resilience [40,41].

4.2. Limitations and Future Work

Several limitations should be acknowledged. First, the use of administrative data introduces potential reporting inconsistencies, particularly among smaller facilities. Second, PSI was constructed using only GDR and NDR, which are closely related mortality measures therefore both indicators may capture overlapping information. Consequently, the significant clustering observed in PSI may largely reflect mortality patterns rather than broader dimensions of patient safety. Third, mortality rates are influenced by factors beyond patient safety, including case severity, referral patterns, hospital type, availability of specialist, and socioeconomic characteristics of the population. Therefore, district with higher mortality rates do not necessarily provide lower-quality care. Fourth, while spatial analysis identifies clustering, it does not explain the mechanisms underlying these patterns. Spatial regression analyses could be conducted to identify factors associated with the observed district-level clustering using variables such as healthcare workforce density, hospital capacity, urban-rural characteristics may help explain why certain districts experience similar PSI outcomes.

Future research incorporating additional patient safety indicators, such as healthcare-associated infection, patient safety incident reports, or complication rates, hospital level and patient-level data would help to address these limitations would provide more comprehensive assessment of patient safety beyond mortality alone. Additionally, future research could consider employing advanced statistical techniques such as regression analysis to understand the underlying causes of high or low hospital performance. Moreover, operational performance and patient safety trends should be investigated over a longitudinal time series to better understand the improvement or decline of regional healthcare. Finally, multiyear assessments should be combined with machine learning models in order to make future-oriented predictions and forecasts of hospital performance. Such an analysis would enable policymakers to foresee future problems and tackle them in a timely manner.

5. Conclusions

This study makes a methodological contribution by implementing comprehensive spatial mapping and spatial autocorrelation analysis within the context of hospital performance and patient safety outcomes. By distinguishing operational performance and patient safety outcomes, the analysis provides a more comprehensive understanding of healthcare quality than single-indicator approaches.

The findings suggest that variations in hospital performance and safety outcomes in East Java are driven primarily by institutional and governance factors rather than geographic proximity. While no significant clustering was observed at the individual hospital level, meaningful spatial patterns emerged at the district level, highlighting the importance of local governance environments, resource distribution, and health system organization in shaping the outcomes.

These results contribute to the growing application of GIS in public health research by demonstrating how routinely collected administrative data can be leveraged to identify geographic disparities and support evidence-informed decision-making. For policymakers, the findings highlight the need to prioritize districts with persistently high mortality while learning from districts demonstrating favorable outcomes. Improving patient safety will require a combination of geographically targeted interventions in high-risk districts and institution-specific strategies tailored to local hospital contexts and capacities.