Study area

Kaohsiung City, located in tropical regions, is Taiwan's third most populous metropolitan area. The tropical monsoon climate, characterized by high humidity and temperature, creates ideal conditions for transmitting the dengue virus [21, 28]. With an average monthly temperature of 25.4 °C and concentrated rainfall of 1968.2 mm from May to September [31], the area is conducive to the breeding of vector mosquitoes, Aedes aegypti and Aedes albopictus, leading to frequent dengue outbreaks in the past 20 years [38].

The study area is focused on the urban core of Kaohsiung, including the Fongshan district and other major metropolitan areas (Fig. 1). The population density in this area is 9097 people/\(}^\), with an average of 20,174 people/\(}^\) and a standard deviation of 12,398 people/\(}^\) across villages. The village serves as the primary unit of analysis, with 528 villages in 12 districts, spanning a total area of 205 \(}^\), and an average size of 0.4 \(}^\) per village. To facilitate visualization and comparisons of hotspot and edge area distributions in the study area, we have highlighted the five regions (Region N1, N2, C1, C2, and S) using blue circles in Fig. 1. These regions were selected based on considerations of population density and the home range of residents within the city.

Fig. 1

Study area. The subgraph on the top-right shows the location of the study area relative to other countries in East Asia. And the main graph shows the main area in Kaohsiung city, Taiwan and the population density in each village. According to the concentration of population density distribution and geographical location, we have hightlighted the five regions (N1, N2, C1, C2, and S) in blue circles shown on the graph in the study area to discuss the hotspot and edge area distribution

Disease data

This study utilizes daily confirmed dengue data obtained from the Taiwan Centers for Disease Control (Taiwan-CDC), covering a period from 1998–2020. This dataset enables the examination of different outbreak scenarios at varying scales. The data fields for each case include the dates of illness onset and notification, the age of the infected individual, the village of residence, and whether the case is indigenous or imported.

Each confirmed dengue case from Taiwan-CDC is tested through laboratory diagnosis. A confirmed case is defined as: (i) positive for dengue virus isolation; or (ii) positive for dengue virus genomic sequence; or (iii) positive for non-structural protein 1 (NS1) by serum antigen test (iv) four-fold increase of dengue virus-specific IgM or IgG antibody in paired serum samples [30].

In this study, we focus on the dengue cases from severe dengue epidemics to examine spatial variations in disease risk. Using data from 1998 to 2020, we identified three particularly severe years: 2002, 2014, and 2015 (represented by red bars in Fig. 2). To examine temporal variations within each epidemic, we divided each year into four periods: the initial period (A), diffusion period (B), peak period (C), and under-control period (D) according to the growth curve of large-scale outbreaks. Each period lasted four weeks, and the number of indigenous cases in each period were aggregated to determine hotspot and boundary detection. The selection of a four-week time period is based on the average generation time of dengue epidemics, typically spans from 16 to 34 days [15], allowing us to capture the disease dynamics and progression effectively. We also defined the epidemic stages as the transition between consecutive periods, including the exponentially-rising stage (A to B), continuously-rising stage (B to C), and declining stage (C to D) (as shown in Fig. 3).

Fig. 2

Yearly dengue indigenous case number in the study area

Fig. 3

An example of dividing weekly data into periods

MethodsStudy framework

The research framework, illustrated in Fig. 4, consists of three key steps. First, the weekly epidemic data from the selected years is divided into four periods. Next, the hotspots and edge areas are detected using a conditional autoregressive (CAR) model. Finally, we assess whether the risk of dengue transmission is higher at the edge areas of a cluster compared to the cluster itself during different stages of each selected year. The difference-in-difference (DID) effect is estimated by comparing the incidence rates between hotspots and edge areas with a spatial error random effect panel model.

Fig. 4Statistical methodsCAR model for hotspot detection

The hotspots are detected using the CAR model [19]. The model utilizes the number of indigenous confirmed cases of dengue in each village, represented as \(_\), to estimate the probability of zero-based point mass distribution (\(_\)) and mean of the mixed distribution (\(_\)). The data likelihood model used is a zero-inflated Poisson model, which combines the excess of zero infections in many villages and the count data of dengue cases. It is assumed to be a joint point mass distribution based on the logistic distribution and a Poisson distribution (Eq. 1).

$$_ \sim \mathrm(_,_)$$

(1)

Then, a linear equation, Eq. 2, is used to estimate the number of cases in each village. \(_\) is the intercept, \(\beta\) is the regression coefficient, and \(_\) is the population in \(k\)th village because there would be more cases in a place with more people.\(_\) is given as CAR prior, which stands for a spatial random effect. The CAR prior (Eq. 3) includes only a single set of random effects \(_\). It is given the condition that corresponds to a normal distribution, and a village k is affected by all the neighbors except itself (\(_\)). Here, movement of daily life within short distances is assumed to build the adjacent matrix W. The definition of "neighbor" relies on Rook contiguity, which is underpinned by the assumption that individuals typically move within their neighboring villages in their daily routines. If two villages share the same edge, then \(_=1\); otherwise, if two villages are not adjacent, then \(_=0\). Also, there are two parameters \(\rho\) and \(^\) to be estimated. \(\rho\) represents the parameter of spatial dependence. If \(\rho\)=0, the value between each village is independent; if \(\rho\)=1, there is spatial autocorrelation, and it corresponds to the intrinsic CAR model. \(^\) is proportional to the concentration of this normal distribution.

$$\mathrm\left(_\right)= _\beta + _+ _$$

(2)

$$_| _, W, ^, \rho \sim N(\frac_^__}_^_ + 1 - \rho }, \frac^}_^_ + 1 - \rho })$$

$$^\sim \mathrm-\mathrm(1, 0.01)$$

$$\rho \sim \mathrm(0, 1)$$

Hotspot probabilities (\(}_\)) are estimated straightly by \(_}^\), meaning the probabilities, which are estimated cases in each village, are greater than \(_\) (Eq. 4). And \(_\) twice the average cases in each village is used as the threshold [35] in the calculated period. From the posterior distribution p(\(_\left|\mathrm\right)\), there would be G samples \(_}^\) (g = 1…G), indicating the predicted case number of each village i after computing the CAR model by a Monte Carlo Makov Chain (MCMC). In this study, G is set up to be 100, so each village has 100 estimated case numbers. Next, a village with a hotspot probability greater than or equal to P (P≧0.7) is defined as a hotspot area, and its neighbors are called the edge areas.

$$}_\equiv \widehat\left(_>_|y\right)= \frac_}^>_}$$

(4)

Bayesian areal Wombling for boundary detection

The boundary detection method was first introduced by Womble [37] to measure the rate of species movement. Later, it was expanded upon by Barbujani et al. [2] and Bocquet-Appel and Bacro [5] for point or grid data analysis. For taking spatial autocorrelation and uncertainty into account, a Bayesian areal Wombling approach was developed using an intrinsic conditional autoregressive (CAR) prior as a spatial random effect [7, 22, 23, 33].

Bayesian areal Wombling is used to detect the boundary on the edges of the hotspot areas, also based on the CAR model as mentioned previously. From the posterior distribution of the CAR model, the absolute difference value between each village of the hotspot area and its each edge area neighbor would be the boundary likelihood value (\(_\)), which is the case difference between two neighboring villages between a hotspot area and an edge area (Eq. 5).

$$_}^=\left|_}^-_}^\right|$$

(5)

$$}_\equiv \widehat\left(_>_|y\right)= \frac_}^>_}$$

(6)

Boundary probabilities (\(}_\)) are estimated by counting the number of \(_}^\) in G samples that exceed a threshold \(_\), meaning that the difference between each pair of neighbors is more than \(_\) cases (Eq. 6). For boundary detection, twice the average cases in each village is also used as the threshold in the boundary probability. Each edge between a hotspot area and an edge area would obtain a boundary probability. Then, an edge area village that shares a boundary probability greater than or equal to P (P≧0.7) with a hotspot is defined as a significant edge area (Fig. 5). The remaining edge areas would be defined as non-significant edge areas if they are not categorized as significant edge areas. The distinction between significant edge areas and non-significant edge areas lies in their proximity to hotspots. Significant edge areas exhibit a statistically significant disease boundary with neighboring hotspots, resulting in a notable difference in confirmed cases. This aims to observe whether there is a difference between each edge area adjacent to hotspots, and two types of edge areas.

Fig. 5

Definition of the 3 areas. The left is a schematic graph of incidence choropleth map, the darker the red it, the more confirmed cases, and vice versa. The right is a schematic graph for delimiting hotspot area, significant edge area, and non-significant edge area

Difference-in-difference estimation (DID)

The difference-in-difference (DID) method is a commonly used technique in quasi-experimental research to assess the impact of specific changes, such as a new policies or medical treatments, across multiple groups and time periods [1, 36]. This design involves two groups (control and treatment) and two time periods, assuming that the trends of the groups would be similar without any intervention or treatment [10, 36]. In our study, we make the assumption that the trends in different areas are the same. The underlying rationale behind the "same trend" assumption is to designate one area, particularly the hotspot area, as a reference point for comparison. This approach enables us to assess whether other areas exhibit a trend of increasing rates exceeding that of the hotspot area, thereby identifying a significant DID effect.

DID analysis involve two differences: the difference in time or the natural trend of the incidence rate from one period (\(}_\)) to another (\(}_\)) (as shown by the difference between the black dotted lines in Fig. 6), and the difference in the incidence rate between the two areas (constant difference in Fig. 6). The final result of the DID analysis is the difference between the original trend in the period \(}_\) and the specific area (difference in difference in Fig. 6). Figure 6 demonstrates the situation in which the disease trend is either increasing or decreasing, with a positive value of the DID estimator.

Fig. 6

Concept of difference-in-difference

Following the identification of hotspot areas, significant edge areas, and non-significant edge areas in each period, we applied DID regression to compare the two area groups. It is to determine if the incidence rate in the edge area is significantly higher than in the hotspot area and to compare differences between the two types of edge areas. The DID regression utilizes the spatial error random effect panel model, as described in the following equations (Eqs. 7–9).

$$y=_+(_}\otimes GROUP)_+(TIME\otimes GROUP)_+u$$

(7)

$$u= \rho (}_}\otimes _)u+\varepsilon$$

(8)

$$\varepsilon =(_}\otimes }_)\mu +v$$

(9)

The spatial error random effect panel model proposed by Kapoor et al. [17] considers potential spatial autocorrelation in spatial data. The dependent variable, y, is a vector of incidence rates in each village at two time periods (Eq. 7). GROUP is a dummy variable comparing two of the three area groups, and TIME is a dummy variable indicating the pre- or post-period. The coefficients \(_\), \(_\), and \(_\) represent the intercept, the total difference between area groups, and difference-in-difference between time and area. \(}_}\) is an identity matrix of dimension T, and \(_\) is a spatial weights matrix in an N × N form built by Rook contiguity. The disturbance vector, u, contains two terms—a spatial autoregressive parameter (ρ) ranging from -1 to 1, and cross-sectional specific effects (\(\mu\)). \(_}\) is a T × 1 vector of ones, and \(}_\) is a N × N identity matrix, and \(v\) is the innovation varying over different villages and periods. \(_\), \(_\) and \(_\) should all be independent and identically distributed (i.i.d). The main result of concern is the significance of \(_\) (DID estimator), which indicates if the change in incidence rate in one area is greater than another.

The DID estimator (coefficient \(_\)) measures the growth rate of risk by comparing the incidence rate between two groups (hotspot and significant edge area) over two time periods (\(}_\) and \(}_\)). The model determines if the change in incidence rate in one area is statistically greater than in another. The comparison is made between the incidence rates in the same areas in both periods, with the significant edge area represented as 1 (GROUP = 1) and the hotspot area as 0 (GROUP = 0). The main result of interest is the significance of the DID estimator.

All statistical analyses were performed in R version 4.2.1. The “CARBayes” package was used for CAR models and Bayesian areal Wombling [18], and the “splm” package was used for DID estimation [25].

Sensitivity analysis of parameter settings

We conducted sensitivity tests for identifying hotspots and boundaries with different thresholds and probabilities. First, the values for the boundary and hotspot detection thresholds (c) were set in the range of 1.0 to 3.0, with a fixed probability of boundary (\(}_\)) and hotspot (\(}_\)) of 0.7 (P = 0.7). These thresholds correspond to a multiple (c times) of the average cases in the study area during a specific period. Second, we varied P from 0.5 to 0.9 with fixed thresholds (c) at 2.0 to determine the robustness of results under different parameter settings.