Exploring Dengue Dynamics: Multi-Scale Analysis of Spatio-Temporal Trends in Ibagué, Colombia

1. Introduction

Dengue fever (DENV), a mosquito-borne viral disease [1], has become a critical public health issue globally, particularly in tropical and subtropical regions [2,3]. Ibagué, Colombia, exemplifies an urban area significantly affected by dengue, with its incidence rising notably over the past years [4]. Ibagué's rapid urbanization over the past two decades has led to densely populated, low-income neighborhoods that often lack regular access to water and adequate infrastructure. These socioeconomic conditions have been linked to higher dengue incidences and mosquito populations [5,6,7]. Additionally, the city's specific environmental characteristics, including its elevation (1,225 meters above sea level) and average temperature (24°C), create a conducive environment for Aedes mosquitoes, further exacerbating the situation [8,9,10].

Dengue transmission dynamics are complex and influenced by a myriad of factors, including urbanization, demographic changes, and environmental conditions [11,12]. Prior research has demonstrated the important role of socioeconomic and demographic variables in dengue spread [13,14]. However, a significant gap still needs to be addressed in understanding the interaction of these variables at different urban scales, namely levels of aggregation, particularly in rapidly urbanizing cities in developing countries [15].

Different approaches have been considered to identify and understand the underlying behavior of diseases in urban environments and their relationship with socioeconomic and environmental variables. For instance, in 2001, authors studied the relationship between neighboring socioeconomic effects and health outcomes, finding consistent results on the influence of neighborhood effects on health [16]. Other studies between 2000 and 2010 asserted a correlation between neighborhood distribution and health status, suggesting not only the use of distance and proximity to determine neighborhoods’ composition but also random effects regarding the neighboring structures to account for the possible connectivity between them [17,18,19,20]. The latter has been widely used in the statistical modeling of DENV using a spatially structured random effect to consider spatial correlation [21,22,23].

To assess the significance of socioeconomic and environmental variables regarding dengue incidence and relative risk (RR), authors have implemented a General Linear Model (GLM) Log-linear model to associate the socioeconomic typology with the risk of dengue infections during winter in a 250 x 250 m grid of Delhi, India [24]. They found out that the lack of constant access to tap water was a risk factor for dengue infections; also, densely populated areas did not necessarily have higher mosquito populations, mainly because of the need for available breeding sites. Finally, high DENV seropositivity was found in wealthier neighborhoods, while low mosquito exposure was explained by mobility within the city. However, this study only considered socioeconomic variables but not demographic predictors, and they ignored possible random effects regarding spatial and seasonal variability. Other authors followed a different approach, implementing a Bayesian Hierarchical Model to determine the association between the relative risk of dengue and lag response in hydrometeorological hazards at a microregion level. The results showed that droughts increase the relative risk of dengue infection in urban areas, while wet conditions increase the risk of infection in rural areas [25].

Specifically in Colombia, authors have implemented a similar methodology for DENV in the city of Cali, using a Space-Time Conditional Autoregressive Model with added autocorrelated random effects for spatial structures and time [26]. This study was conducted on a neighborhood level, where socioeconomic variables were controlled using a Principal Component Analysis (PCA) approach, and several environmental variables were included and lagged to consider the delayed effect. The results showed that lagged weather variables could help to identify when the peaks in the risk of transmission occur. Additionally, they proved that dengue infections are not exclusive to poor areas, and the risk of infection is related to spatial and temporal distribution. The proposed aggregation level of neighborhoods offers very sparse data observations with clear socioeconomic and demographic trends.

Bayesian models have shown that they can be used to assess significance and make inferences about the predictors [22,25,26,27]. The implementation usually calculates the posterior distributions via Markov Chain Monte Carlo Simulations (MCMC) or Integrated Nested Laplace Approximation (INLA). MCMC has been widely used to identify marginal distributions and inference [28,29]; however, it takes a very high computational time [30]. To solve this issue, authors proposed the INLA methodology, which uses local approximations and Laplace transformations, providing very similar results around 300 times faster.

While the importance of spatial and temporal variables in dengue transmission is recognized [31], limited research has been conducted on integrating these factors at different levels of urban spatial aggregation. This study aims to bridge this gap by leveraging detailed demographic and socioeconomic data from the census, provided by the National Administrative Department of Statistics (Departamento Administrativo Nacional de Estadística - DANE), and environmental variables derived from satellite imagery and previous studies. We seek to unravel the spatial and temporal dynamics influencing dengue incidence in Ibagué between 2013 and 2018, examining these factors across various levels of spatial aggregation - Manzanas, Secciones, Sectores, and Comunas [32]. We introduce a novel approach by employing Geographically Weighted Regression (GWR) to isolate key socio-economic and demographic predictors at varying spatial scales [33]. Additionally, the use of INLA models allows for an in-depth examination of spatio-temporal correlations and their posterior distributions [34], offering new insights into the localized dynamics of dengue transmission.

The methodological innovation of this study lies in its tripartite modeling strategy, encompassing spatial, temporal, and combined spatio-temporal models at each level of spatial aggregation. This approach allows for a comprehensive analysis of the varying impacts of different factors on dengue incidence, providing a nuanced understanding of the disease's transmission dynamics in an urban setting. Our findings will contribute significantly to public health, particularly in developing targeted dengue control strategies [35]. The study's framework also offers a valuable model for similar epidemiological investigations in other urban settings, enhancing our understanding of vector-borne diseases in global urbanization trends.

2. Materials and Methods

2.1. Study Site and Data

Ibagué is the most populated city in the department of Tolima, with an estimated population of 541,101 people for 2020 [36]. According to DANE, the city’s urban area is divided into 13 Comunas, 65 Sectores, 272 Secciones, and 5,574 Manzanas, which will be used as levels of spatial aggregation. Groups of the minimum level, Manzanas, create the rest of the levels, as can be observed in supplementary Figure S1. The demographic and socioeconomic data was obtained in all four levels from the National Geostatistical Framework (Marco Geoestadístico Nacional - MGN) and National Population and Dwelling Census (Censo Nacional de Población y Vivienda - CNPV), which were compiled in 2018 [37].

Socioeconomic predictors included variables related to strata, water access and disposal, access to gas, garbage pick-up service, and internet connection, which allows to identify lower income and utility access areas inside the city. Demographic variables divide the population according to age, gender, and educational level. Finally, environmental variables were obtained from satellite images via Google Earth Engine using MODIS1A1 for Mean Temperature [38], MODIS13Q1 for NDVI [39], CHIRPS PENTAD for Total Precipitation [40], and previous studies on the city for wet days and hot days (over 32°C) [41]. The variables and their description are shown in supplementary materials.

Data on dengue cases was obtained from the dataset provided by the local government. Out of the total 17,707 DENV recorded from 2013 to 2018, 16,183 were included according to the spatial levels of the city. The remaining cases had no readable address or were reported in the city's rural area. Figure 1 displays the evolution of the virus over the years. 2013 and 2015 showed the most extensive outbreaks, with 5,383 and 4,885 cases respectively, with another outbreak following during the first months of 2016. Notably, the number of cases in 2017 and 2018 was significantly lower, with each year reporting fewer than 1,000 cases.

2.2. Methodology

2.2.1. Variable Selection and Transformation

To demonstrate association between the variables and dengue incidence per 100,000 inhabitants, and reduce the dimensionality of spatial variables, we used GWR through the package GWModel in R [42,43]. This methodology allows the creation of different local regressions using Ordinary Least Squares (OLS) for every spatial feature, adding a weight parameter that was obtained from the distance between the geometries and a calculated bandwidth [44,45]. This analysis was performed for socioeconomic and demographic predictors on every level of aggregation. Only variables that were shown to be significant, at least in one global regression at any level of aggregation, were used as final predictors for the inference model.

We also used a Wavelet Coherence analysis to confirm an existing correlation between DENV cases and temporal variables. This was used to analyze non-stationary time series. The methodology implements a decomposition between time and frequency using a windowed Fourier transform, which allows for local time-frequency properties while adjusting for high and low-frequency structures [46]. The wavelet coherence was computed for the aggregated environmental variables in the whole city. Results are evaluated graphically using the biwavelet library in R, keeping the variables that showed high correlation and significance in the wavelet coherence analysis as the final predictors [47]. All temporal predicted were included, since they displayed high correlation with DENV cases as shown in Figure S2.

The temporal predictors were later lagged using a Distributed Lag Nonlinear model (DLNM). This methodology considers delayed effects and nonlinear relationships between dengue and environmental time-dependent variables [48]. The package dlnm was used on R [49]; obtaining, as a result, a matrix that accounts for the nonlinear exposure and a delayed effect.

2.2.2. Model fitting.

The latent marginal distribution of chosen predictors was approximated using INLA, wherein the spatio-temporal influence of these predictors was considered. This approach relies on Latent Gaussian Models (LGM), specifically, a Latent Gaussian Markov Random Field with a sparse and factorizable precision matrix. Such a structure enables numeric approximations, leading to quicker outcomes than conventional LGM techniques like Markov Chain Monte Carlo (MCMC) [34].

In this project, we specifically fitted a negative binomial model to the number of cases $Y_{st}$ for each spatial unit $s$ at a given time $t$, having an estimated mean of ${\mu }_{st}$ and a dispersion parameter $\varphi$. This model allows us to account for overdispersion in the number of cases at each scale. Our link function considers ${\mu }_{st}$ from the population $p_{st}$ and the monthly incidence at the same location and time $\left({\rho }_{st}\right)$. The model is depicted in equations (1) to (3).

$$Y_{st}|{\mu }_{st}~NB({\mu }_{st},\varphi )$$

(1)

$$\log \left({\mu }_{st}\right)=\log \left(p_{st}\right)+\mathrm{l}\mathrm{o}\mathrm{g}\left({\rho }_{st}\right)$$

(2)

$$\log \left({\rho }_{st}\right)=X^T\beta +{\gamma }_{st}+{\eta }_{st}$$

(3)

We also included two random effects that were considered to account for unobserved variability. An unstructured random effect for seasonal autocorrelation, regarding possible relationships in time for each structure along the months, being cyclic over the six-year analysis; and a second structured random effect, that encompasses spatial autocorrelation between the neighborhoods along the years regarding interconnection, interventions, herd immunity, etc. [25] The incidence was calculated from the fixed effects $X^T$ and the two random effects: ${\gamma }_{st}$ for the unstructured effect using a random walk, and ${\eta }_{st}$ for the structured one using a Besag-York-Molliè model [50]. Precision priors are defined from the precision parameter $P\left(\sigma >0.5\right)=0.01$.

An adaptative strategy was selected as it is considered the best fit from Gaussian and simplified Laplace approximations. Finally, the hyperparameters posterior distributions were calculated with a Central Composite Design as it offers the best tradeoff between precision and computational time among the possible strategies implemented in the R library INLA [30,34,51].

Three models were adjusted per level: one containing only spatial variables (socioeconomic and demographic), a second containing only temporal variables (environmental), and a third containing both. This enabled us to distinguish between the effects of spatial and temporal covariates, while assessing whether incorporating both yielded a more informative model. Comparison was performed for the models at each level using the Deviance Information Criterion (DIC) and Widely Applicable Information Criterion (WAIC).

2.2.3. Aggregation Level Comparison

From selected best-fitting models at each level of aggregation, fitted values were obtained from the marginal posterior distributions of the selected models at each level of aggregation. These values were then compared with the observed data over the study period to determine the best fitting model at each level considering the Root Mean Squared Error (RMSE).

3. Results

3.1. Variable Selection

The Geographically Weighted Regression (GWR) analysis revealed a significant association for seven critical variables at a minimum of one spatial level over different years. These variables encompassed critical aspects of the socio-economic landscape, including population density (Density), sewage connection (Sewage), gas connection (Gas), garbage collection service (Garbage), population with higher education degrees (Higher Ed.), percentage of women (Women), and percentage of children (Children).

Figure 2 visually presents the space-time trends of these variables across diverse spatial levels. The visualization aids in understanding the spatial patterns and temporal changes exhibited by each variable, contributing to a comprehensive understanding of the intricate dynamics within the studied region. This provides a concise summary of the outcomes derived from the GWR analysis, emphasizing the spatial and temporal significance of the identified variables.

3.2. Model Fitting and Comparison

Table 1 shows the results of model selection based on DIC and WAIC by level of aggregation. It is important to note that these models can be compared within the same level but not across different levels. A significant shift in both DIC and WAIC values is evident in the temporal models as opposed to the spatial models. Notably, the spatio-temporal models, which integrate a comprehensive set of covariates encompassing socioeconomic, demographic, and transformed environmental predictors, exhibit the most favorable scores in these comparison metrics across all levels and were selected as the best model for multilevel comparison.

Figure 3 provides an insight into the spatial fixed effects of each spatiotemporal model inferred from posterior distributions, showing the variations of these effects across different levels of spatial aggregation.

At the level of Comunas, spatial variables exhibit non-significance, as indicated by the inclusion of zero within the 95% credible intervals for all variables. This lack of significance aligns with the comparable DIC and WAIC values observed for both temporal and spatiotemporal models in Table 1, implying that incorporating spatial variables at this aggregation level does not substantially enhance the model’s explanatory capacity for DENV incidence.

For the intermediate level of Sectores, certain variables such as garbage collection, higher education, and the percentage of children demonstrate significance and display inverse correlations with DENV incidence. The negative coefficients suggest that improved access to garbage collection, higher educational attainment, and a larger proportion of children are associated with reduced DENV spread, potentially highlighting the impact of enhanced public services and education on disease mitigation.

At the Secciones level, most spatial variables exhibit significance, with garbage collection and the percentage of women showing negative correlations with DENV cases. This inverse relationship suggests that areas with more efficient waste management and a higher proportion of women tend to have lower DENV prevalence. Moreover, factors such as higher population density, gas connections, education levels, and the percentage of children consistently demonstrate negative correlations with DENV cases, echoing the trends observed at the Sectores level and underscoring the influence of these variables on disease incidence.

Finally, at the granular level of Manzanas, sewage connection emerges with a unique positive correlation with DENV cases. Conversely, population density, higher education, and the percentage of children maintain inverse correlations with DENV incidence, aligning with observations made at the Secciones level. This consistency across different levels of spatial analysis suggests that certain factors consistently relate to lower disease incidence despite the finer granularity of data.

Integrating this insight into the preceding analysis highlights a consistent negative correlation between spatial variables such as higher education, population density, and the percentage of children with DENV incidence across various levels, while the association between garbage collection services and DENV cases appears less definitive, suggesting disparities in waste management service.

Complementing the spatial analysis, Figure 4 introduces lagged temporal predictors through contour plots, which offer insight into the temporal dynamics of the disease. The data underlying the cyclical influence of weather patterns on the RR associated with DENV. For instance, temperature-related metrics, such as mean temperature and the number of days exceeding 32°C, reveal a lower RR at cooler temperatures and greater lags, transitioning to a higher RR as temperatures climb and lag decreases. Conversely, precipitation indicators, such as total precipitation and number of wet days, demonstrate an inverse relationship, with higher precipitation levels correlating with a decreased RR in subsequent periods. Additionally, the NDVI exhibits variability and lacks consistency across different levels of aggregation, suggesting complex interactions between vegetation density and disease transmission that warrant further investigation.

3.3. Level Comparison

Figure 5 depicts the compared values of observed DENV cases and the fitted values from posterior distributions at each level, aggregated for the entire city.

The model's performance varies across different levels of spatial aggregation. At the Comunas level, the model achieves a high degree of correlation with the actual DENV case data, although it tends to underestimate case numbers during outbreak peaks, a trend that was particularly pronounced towards the end of 2015. In times of lower disease incidence, such as December 2017 and June 2018, the model also falls short of accurately capturing the case numbers. Despite these limitations, it performs commendably in periods of low case counts, notably, throughout most of 2017 and the early months of 2018.

Moving to the Sectores level, the model has a propensity to overestimate the incidence of DENV cases, with this trend being especially evident in January 2015 and January 2016. This tendency for overestimation continues through the endemic years of 2017 and 2018. Moreover, there is a noticeable misalignment in the timing of predicted outbreaks compared to actual data, highlighting a phase discrepancy between model predictions, and observed case trends.

Similarly, in Secciones level, there is a need for phase alignment. Still, the model demonstrates an accurate fit during epidemic periods, both in terms of case count and pattern, as seen between September 2015 and January 2016. This level also accurately captures smaller peaks during endemic periods, such as December 2016 and June 2018. However, during periods characterized by generally low DENV activity, such as the late 2017 to early 2018 timeframe, the model tends to overpredict the number of cases, indicating a challenge in accurately modeling low incidence rates.

At the most granular level of spatial aggregation, Manzanas, the model's fit shows the greatest fluctuation among all levels. While it aligns more closely with the observed data during epidemic periods, similar to Secciones, its predictive performance is significantly less accurate during endemic years, indicating a disparity in model fit across different disease prevalence periods.

The differences in model performance underscore the challenges in capturing the complex dynamics of DENV transmission, which vary temporally with epidemic and endemic cycles and spatially at different levels of urban granularity. These insights emphasize the need for models that can adjust to both the scale of analysis and the fluctuating nature of disease transmission, highlighting the intricate balance between spatial resolution and predictive accuracy in epidemiological modeling.

Table 2 rectifies the observed behavior observed in Figure 5. Lower RMSE was found for Comunas, since it provided the best overall fit; however, the level of Secciones provided a lower RMSE than Manzanas and Sectores, which might be due to the better fit exhibited during epidemic seasons.

4. Discussion

GWR analysis shows variation for spatial covariates among levels of aggregation. Socioeconomic predictors are mostly significance at the lowest level, Manzanas, while demographic predictors were also significant at intermediate levels like Secciones or Sectores. Only two predictors were significant at Comunas. GWR analysis solely identifies significant spatial variables and overlooks non-linear interaction, which is crucial for understanding the endemic-epidemic patterns of DENV. This limitation may result in models that only partially capture the disease's dynamics.

Our study shows that the spatio-temporal models, which integrate spatial and temporal variability, generally offered a better fit, and were selected for further analysis at all levels of aggregation. However, the improvement in the spatio-temporal models at levels like Comunas and Sectores was minimal, indicating that the added complexity of spatio-temporal models may only sometimes lead to significantly improved fit. This underscores the importance of careful consideration when increasing model complexity and highlights the need to balance detailed spatial-temporal dynamics with model simplicity.

The significance of spatial predictors in DENV incidence varied markedly across different levels of spatial aggregation. At the highest levels, namely Comunas and Sectores, many variables were found to be non-significant, possibly due to the homogeneity within these broader spatial categories. The lack of variability within spatial covariates at these levels leads to a limited ability to discern significant impacts on DENV cases, as observed in the narrow covariate ranges detailed in supplementary Table S1. In contrast, population density was notably significant at the more granular Secciones and Manzanas levels and inversely correlated with DENV incidence. However, This counterintuitive finding is supported by previous studies that have suggested that higher densities may not favor mosquito breeding, particularly in areas with sufficient sanitation and utility services [52,53].

Socioeconomic variables, such as sewage, gas connection, and garbage collection service, exhibit diverse correlations with DENV cases. Sewage connection is positively correlated, attributed to urban infrastructure and vector ecology (adaptation to breed in manmade environments), including improperly designed systems creating mosquito breeding sites [54,55,56,57]. Conversely, areas with higher gas connection rates tend to have lower DENV incidence, reflecting socioeconomic status. The relationship with garbage collection services varies; a positive correlation at the Secciones level and an inverse relationship at Sectores suggest complex dynamics influenced by local practices and infrastructure. [58,59].

The observed variability among socioeconomic variables underscores the intricate interplay of individual and collective dynamics, sometimes resulting in counterintuitive outcomes. Nevertheless, these variables offer valuable insights into the internal dynamics of spatial distribution, exhibiting distinct characteristics across different levels of spatial aggregation.

Demographic predictors, including higher education, percentage of children, and percentage of women, demonstrate a consistent pattern across spatial levels. Higher educational attainment in populations may lead to increased implementation of disease prevention measures, potentially reducing breeding sites and subsequent DENV cases [60,61]. The inverse correlation with the percentage of children may reflect the demographic profile of DENV cases during specific outbreaks, with adults and young adults being more frequently affected. Lastly, the percentage of women may reflect the broader demographic composition of the city and the roles women play in household management and potential exposure to mosquito breeding sites.

Variations in correlation between spatial variables and DENV incidence across different spatial aggregation levels indicate diverse roles of socioeconomic and demographic factors in disease transmission dynamics. While broader scales may obscure these factors amidst other influences, finer scales highlight their heterogeneity, allowing for a more precise understanding of their impact on DENV transmission. This underscores the complexity of disease incidence modeling and emphasizes the importance of considering scale when interpreting variable influences.

Despite consistent trends in temporal variables across all spatial aggregation levels, subtle variations in relative risk suggest differing degrees of correlation with DENV cases. Temperature and precipitation influence were relatively consistent across the city's spatial structures, resulting in similar impacts on relative risk regardless of aggregation level. Our results show an association between Mean Temperature and increased Relative Risk (RR) during high-temperature periods (zero-month lag), and elevated disease risk due to total precipitation at a three to six-month lag across Comunas, Sectores, and Secciones, which underscores the influence of weather patterns on DENV transmission in urban areas. This relationship, consistent with regional climate patterns where cooler seasons precede warmer periods by approximately four to six months [62], which is further intensified by ENSO phenomena, such as the El Niño events of 2015 and 2016. These extreme weather events, characterized by higher temperatures and drought, impact mosquito breeding and access to utilities in vulnerable areas, leading to community adjustments in water management practices and influencing DENV transmission dynamics [25,63].

Days Over 32°C exhibited a comparable trend with mean temperature, where a higher frequency of hot days at shorter lags was linked to an increased RR. The same relationship is displayed by the number of wet days and total precipitation, with increased RR at two and six-month lags, reflecting bimodal rainy seasons. Finally, NDVI showed irregular patterns, with minimal variation at broader levels like Comunas and Sectores due to its non-seasonal nature.

While the temporal patterns hold consistently across spatial aggregation levels, the changes in RR's for each covariate vary, suggesting differing strengths of correlation with DENV incidence. This variance could stem from how DENV cases are distributed across each spatial level and how spatial covariates account for the observed effects at more granular levels.

This study highlights the importance of considering temporal and spatial variables in understanding DENV transmission dynamics. While temporal variables play a significant role, spatial covariates at finer levels of aggregation are also crucial for a nuanced understanding of DENV transmission. The analysis suggests that the Comunas level model provides the best overall fit for the city, with the Secciones level model closely following. However, the Manzanas level model performs weaker due to extreme case dispersion. Despite this, intermediate aggregation levels like Secciones reveal discernible links with socioeconomic and demographic variables, aiding in understanding local patterns for targeted interventions. This underscores the importance of considering macro and micro-level factors in epidemiological modeling and intervention planning to tailor public health strategies and reduce disease prevalence effectively.

The accuracy of DENV case reporting relies on local population engagement, often hampered by underreporting due to symptom recognition without seeking formal diagnosis and the prevalence of asymptomatic cases [64]. This underreporting significantly impacts case counts' accuracy, hindering model precision [65]. Census data limitations are apparent, with only 2018 data available, assuming socio-economic and demographic variables remain unchanged over six years, neglecting potential variations. Modeling efforts focusing on endemic or epidemic periods may offer immediate insights into socio-economic and demographic influences on disease patterns but may overlook long-term effects. Detailed research at lower observational levels is needed to address data scarcity and the influence of local entomological and virological factors on disease dynamics.