Tuberculosis spatial clustering and district-level determinants in Eastern Province, Zambia: a Bayesian spatial analysis of program surveillance data (2020-2024)

¹Department of Public Health, University of Lusaka, Lusaka, Zambia, ²Medical Department, AIDS Healthcare Foundation, Lusaka, Zambia, ³Department of Pharmacy, Marybegg Health Services, Ndola, Zambia, ⁴Department of Internal Medicine, University Teaching Hospital, Lusaka, Zambia, ⁵Public Health, National Tuberculosis and Leprosy Programme, Lusaka, Zambia, ⁶Department of Public Health, Eastern Province Provincial Health Office, Chipata, Zambia, ⁷Department of Public Health, Copperbelt University, Ndola, Zambia, ⁸Department of Public Health, National Health Research and Training Institute (Formerly TDRC), Ndola, Zambia, ⁹Department of Public Health, Mulungushi University, Livingstone, Zambia

^&Corresponding author
Claude Mumba, Department of Public Health, University of Lusaka, Lusaka, Zambia

Introduction    

Tuberculosis (TB), a preventable and curable disease, affected approximately 10 million people and caused over 1.5 million deaths in 2023, surpassing HIV and malaria as the world's leading infectious killer [1]. Approximately one-quarter of the world's population is estimated to be infected with Mycobacterium tuberculosis, placing millions at ongoing risk of progression to active disease [2]. The burden of TB is disproportionately concentrated in low- and middle-income countries, particularly in sub-Saharan Africa (SSA), where poverty, malnutrition, overcrowding, and high HIV prevalence continue to drive transmission and poor outcomes [3].

Sub-Saharan Africa accounts for nearly 25% of global TB cases and approximately 30% of TB-related deaths [1]. Zambia is among the 30 high TB-burden countries globally and ranks 21^st in terms of TB incidence [1]. The World Health Organization (WHO) estimates Zambia's national TB incidence at 283 cases per 100,000 population, translating to approximately 59,000 new cases annually, with treatment coverage estimated at 93% [1]. Despite progress in TB diagnosis and treatment, persistent challenges such as delayed diagnosis, inadequate contact tracing, underreporting, and health system constraints continue to undermine national TB control efforts [3]. Evidence from national analyses has also highlighted substantial gaps along the TB care cascade, particularly in case detection and linkage to care [4].

Eastern Province is one of Zambia's high-burden regions for TB, characterized by a mix of urban and rural districts with significant cross-border movement with Malawi and Mozambique. These dynamics, coupled with high HIV prevalence, limited healthcare access, and seasonal migration, create an environment conducive to sustained TB transmission [5]. Although TB burden remains high in the province, spatially disaggregated evidence describing the geographic distribution of TB and its underlying determinants remains limited. Most TB studies conducted in Zambia have focused on national-level trends, facility-based outcomes, or clinical characteristics, with limited application of spatial epidemiological methods [4]. Such approaches often overlook geographic heterogeneity and localized clustering of TB, which are critical for guiding targeted interventions in resource-constrained settings. Spatial analytical techniques, particularly Bayesian disease mapping, have emerged as powerful tools for identifying high-risk areas and quantifying geographic variation in disease burden while accounting for spatial dependence and uncertainty [6].

When combined with Geographic Information Systems (GIS), these models allow for visualization of disease patterns and support evidence-based prioritization of interventions, surveillance, and resource allocation. Recent studies from other high TB-burden countries have demonstrated the utility of Bayesian spatial models in identifying TB hotspots and informing targeted control strategies [7]. Therefore, this study used routinely collected TB surveillance data from 2020 to 2024 to examine the spatial distribution of TB incidence and its area-level socioeconomic and clinical determinants in Eastern Province, Zambia. By applying hotspot analysis and Bayesian disease mapping, this study aims to provide localized evidence to support targeted TB control strategies, strengthen surveillance, and guide more precise public health interventions in high-burden settings.

Methods     

Study design and setting: this study was a quantitative secondary analysis of routinely collected TB surveillance data from all 15 districts of Eastern Province, Zambia, covering the period from 1 January 2020 to 31 December 2024. Eastern Province is located in the eastern part of Zambia between 30°-34° east longitude and 10°-15° south latitude, covering approximately 51,293 km². It shares international borders with Malawi to the east and Mozambique to the southeast, while Muchinga, Central, and Lusaka Provinces lie to the north, west, and southwest, respectively [8]. The province comprises 15 districts, with Chipata serving as the provincial capital. According to the Zambia Statistics Agency (ZamStats), Eastern Province had an estimated population of 1.9 million in 2020, of whom 50.5% were female [8]. The province's economy is predominantly agrarian, with substantial cross-border trade and seasonal population movement. HIV prevalence among adults aged 15-49 years was estimated at 7.7% in Eastern Province compared with the national average of 12.3% [8].

Study population: the study population included all TB cases notified through the District Health Information System (DHIS2) across the 15 districts of Eastern Province between 2020 and 2024. Both new and relapse TB cases were included, regardless of age or sex. Exclusion criteria removed cases with unvalidated or incomplete district identifiers, and district-years with reporting completeness below 95%. The analysis spans 1^st January 2020 to 31^st December 2024. Missing data were resolved through verification with district officers, and no districts were excluded as all 15 maintained >95% completeness annually. As this was an ecological analysis, the unit of analysis was the district rather than the individual patient.

Sample and sampling strategy: all 15 districts in Eastern Province constituted the study sample. No sampling was performed, as a census of district-level TB notifications was used. For spatial analysis, districts were weighted according to population size to account for differences in population distribution across the province. GIS-based spatial analysis principles were applied, considering district adjacency and spatial proximity.

Data sources and collection: tuberculosis (TB) notification data were extracted from DHIS2, the official national electronic health surveillance platform managed by the Ministry of Health. The dataset included annual district-level counts of new and relapse TB cases. Socioeconomic and clinical covariates were obtained from secondary data sources, including the Zambia Demographic and Health Survey (ZDHS), ZAMPHIA reports, and ZamStats district profiles. These covariates included poverty indicators, HIV prevalence, malnutrition indicators, and ART coverage at the district level. Geospatial data, including district administrative boundary shapefiles and road networks, were obtained from the ZamStats GIS Division and OpenStreetMap. These datasets enabled spatial linkage and visualization of TB distribution across the province. All datasets were aggregated at the district level and cross-validated against Ministry of Health annual reports for consistency.

Data management: data extraction and preprocessing were conducted using combined automated DHIS2 queries with manual retrieval, Microsoft Excel, PostgreSQL, and R version 4.3, within a defined timeframe. Data cleaning involved ordered identification and resolution of missing values, duplicates, and inconsistencies using data quality dashboards and validation rules. Quality control included cross-checking against facility registers, reconciling with Ministry of Health annual reports, and triangulating with ZAMPHIA and ZDHS, while districts with incomplete data underwent sensitivity checks to confirm completeness thresholds were met. TB notification data were merged with socioeconomic and geographic datasets using standardized district codes. Tuberculosis case counts were aggregated to the district level, and annual TB notification totals were computed. Population estimates were used to calculate expected TB counts and served as offset terms in regression models. Quality assurance included triangulation of TB notification totals with Ministry of Health annual TB reports and verification of spatial coordinate accuracy. The final analytic dataset contained complete data for all 15 districts across the five-year study period.

Study variables: the primary outcome variable was TB notification count, defined as the number of new and relapse TB cases reported per district per year. Independent variables included: i) poverty index, defined as the number of households below the national poverty line per district; ii) HIV prevalence (%), defined as the proportion of adults aged 15-49 years living with HIV; iii) malnutrition rate, defined as the number of undernourished individuals per district, is used as a proxy for food insecurity; iv) ART coverage (%), defined as the proportion of people living with HIV receiving antiretroviral therapy; Population size (log-transformed) was included as an offset term in regression models to account for differences in district population sizes. The calendar year was included to adjust for inter-annual variation. Environmental and migration-related factors were considered qualitatively during interpretation due to the lack of routinely available district-level measures.

Statistical analysis: data were analyzed using Stata version 17.0 (StataCorp, College Station, TX, USA), R version 4.3.1, and ArcGIS Pro version 3.0 (Esri, USA). Descriptive statistics were used to summarize TB notifications and district-level socioeconomic and clinical characteristics. Means and standard deviations were calculated for continuous variables. Choropleth maps were generated in ArcGIS Pro to visualize the spatial distribution of TB notifications, poverty, HIV prevalence, and malnutrition across districts, providing an initial spatial overview before inferential analysis. Statistical methods were selected based on the hierarchical nature of the data (districts nested within the province), the count nature of the outcome variable (TB notifications), and the presence of spatial dependence among neighbouring districts. Given the risk of Type I error linked with multiple comparisons in the spatial analysis, we applied a False Discovery Rate (FDR) adjustment to the Getis-Ord Gi* hotspot results (Benjamini-Hochberg procedure). The FDR-adjusted p-values were used to identify statistically significant hotspots, with an FDR threshold of 0.05. For the Poisson regression analysis, we restricted the primary inference to pre-specified covariates based on a priori hypotheses (HIV prevalence, ART coverage, poverty, malnutrition, and urban/rural status), thereby limiting the number of statistical tests and controlling the family-wise error rate.

Poisson regression analysis: associations between district-level TB notifications and selected determinants were assessed using Poisson regression models with a log population offset. Cluster-robust standard errors were applied to account for the potential correlation of observations within districts over time. The model specification was:

log(λ_i) = log(Population_i) + β₀ + β₁(Poverty_i) + β₂(HIV_i) + β₃(Malnutrition_i) + β₄(ART_i)

Where (λ_i) represents the expected number of TB notifications in district i. Incidence rate ratios (IRRs) with 95% confidence intervals (CIs) were reported. Model diagnostics included assessment of multicollinearity using variance inflation factors (VIFs).

Spatial hotspot analysis: spatial clustering of TB notifications was assessed using the Getis-Ord Gi* statistic implemented in ArcGIS Pro. A distance-based spatial weights matrix was applied to define neighboring districts. Statistically significant hotspots were identified by positive Z-scores greater than 1.96, while cold spots were identified by negative Z-scores less than -1.96, at a significance level of p < 0.05. Results were visualized using a red-blue gradient, where deep red represented hotspots and deep blue represented cold spots The Gi* statistic was computed as:

Where: x_j = TB notification count in district j; w_ij = spatial weight between districts i and j; n = total number of districts; X̄ = mean TB notification count across all districts; S = standard deviation of TB notification counts.

Bayesian disease mapping: to account for spatial autocorrelation and estimate district-level TB risk, a Bayesian hierarchical Poisson disease mapping model based on the BYM formulation was fitted using the Integrated Nested Laplace Approximation (INLA) approach implemented in the R-INLA package. Let Y_i denote the observed number of TB cases in district i. Likelihood: Y_i ∼Poisson(λi).

The log of the expected TB count was modeled as:

Log(λi) = log(E_i) + α + β1X_1i + β₂X_2i + β₃X_3i + ui + vi.

Where: E_i = expected number of TB cases in district i, calculated as the district population multiplied by the provincial mean TB notification rate (offset term); α = overall intercept; X_1i = scaled poverty rate; X_2i= HIV prevalence (%); X_3i= scaled malnutrition rate; u_i = spatially structured random effect capturing spatial dependence between neighboring districts (conditional autoregressive component); v_i = unstructured random effect accounting ↓ for district-specific heterogeneity.

Spatial random effects specification

Structured effect (CAR prior)

Where: j∼_i = indicates neighboring districts; n_i = number of neighbors of district i; τ_u = spatial precision parameter.

Unstructured effect

v_i ∼Normal(0, τ_v^-1)

Priors: i) βk ∼Normal(0,10); ii) τu ∼Gamma(1,0.01); iii) τv ∼Gamma(1,0.01).

Relative risk interpretation

RR_i = λ_i / E_i

Districts with: i) RR_i > 1 → higher-than-expected TB risk, ii) RR_i < 1 → lower-than-expected TB risk.

Posterior distributions were estimated through INLA. Relative risk (RR) values with 95% credible intervals (CrI) were computed for each district and visualized as smoothed risk maps. Model convergence and goodness-of-fit were evaluated using the Deviance Information Criterion (DIC), Watanabe-Akaike Information Criterion (WAIC), and posterior predictive checks.

Bayesian manual analysis

For conceptual illustration, a simplified non-spatial Poisson regression model was specified as: likelihood: Y_i ∼Poisson(λ_i)

Log-linear model with population offset:

Log(λ_i) = log(Population_i) + β₀ + β₁x_1i + β₂x_2i + β₃x_3i

Standardization of predictors. All predictors were standardized as:

Where: i) X̄_k = mean of predictor k ii) SD(X_k) = standard deviation of predictor k.

Poisson probability mass function.

Log-posterior for regression coefficients. The log-posterior distribution was expressed as:

Where: β = (β₀, β₁, β₂, β₃); and P(β_j) denotes the prior distribution for coefficient β_j.

Model assessment and sensitivity analysis. Model adequacy and robustness were evaluated using: i) DIC, ii) WAIC, iii) posterior predictive checks.

Sensitivity analyses included: i) re-estimation after excluding districts with extreme TB counts (e.g., Petauke), ii) comparison of BYM results with alternative Poisson-lognormal spatial models, iii) evaluation of parameter stability under alternative prior specifications. Findings were consistent across all model variants, confirming robustness of spatial associations.

Data visualization: spatial outputs were visualized using GIS, including: i) TB notification rates, ii) hotspot and cold-spot clusters (Getis-Ord Gi*), iii) Bayesian smoothed RR maps. All maps were projected using the WGS84 coordinate system to ensure geographic comparability. Model convergence and goodness-of-fit were evaluated using the Deviance Information Criterion (DIC), Watanabe-Akaike Information Criterion (WAIC), and posterior predictive checks. The DIC for the final BYM model was 89.3, and the WAIC was 90.1, indicating adequate model fit relative to the null model (DIC = 112.5).

Ethical considerations: ethical approval was obtained from the University of Lusaka Biomedical Research Ethics Committee (Reference No: FWA00033228-680/2024) and the National Health Research Authority (NHRA/2542/29/07/2025). Permission to use DHIS2 data was granted by the Ministry of Health and the Eastern Province Provincial Health Office. The study used aggregated, de-identified data, and no individual-level identifiers were accessed.

Results     

Descriptive epidemiology of tuberculosis in Eastern Province (2020-2024)

A total of 14,145 TB notifications were reported across the 15 districts of Eastern Province between 2020 and 2024. TB notification counts varied substantially across districts, ranging from 39.8 to 595 cases, with a mean of 188.6 cases per district (SD ± 158.1). The mean district population was 179,794.5 (SD ± 72,576.8), with population sizes ranging from 55,504 to 355,306. The mean poverty index score was 4,441.2 households per district (SD ± 2,273.4), while the mean HIV prevalence rate was 4.25% (SD ± 2.38%), ranging from 1.8% to 9.5%. The mean ART coverage among people living with HIV was 88.47% (SD ± 3.64%), with values ranging from 81.24% to 97.47%. The malnutrition rate exhibited considerable variability across districts, ranging from 137 to 2,922 individuals, with a mean of 720.3 (SD ± 714.8), as shown in Table 1.

Determinants of TB incidence (poisson regression)

Results from the Poisson regression analysis with a population offset and cluster-robust standard errors are presented in Table 2. Urban districts recorded significantly higher TB notification rates compared to rural districts (IRR = 1.73, 95% CI: 1.53-1.95, p < 0.001). HIV prevalence was positively associated with TB notifications, with a one-unit increase corresponding to a 4.5% increase in TB notification rates (IRR = 1.045, 95% CI: 1.022-1.069, p < 0.001). ART coverage was also positively associated with TB notification rates, with a one-unit increase associated with a 13.7% higher notification rate (IRR = 1.137, 95% CI: 1.122-1.152, p < 0.001). This substantial effect likely reflects enhanced TB screening and case finding in high-ART-coverage settings, rather than a causal effect of ART increasing TB incidence. The poverty index showed a statistically significant association with TB notification rates; however, the magnitude of the effect was small (IRR = 1.00014, 95% CI: 1.000111-1.000173, p < 0.001). The malnutrition rate was not statistically associated with TB notification rates in the Poisson regression model (IRR = 0.99999, 95% CI: 0.99992-1.00005, p = 0.686) as shown in Table 2. Further, a negative binomial regression model was fitted as a sensitivity analysis to assess the robustness of findings under potential over-dispersion in the data. The results are shown in Annex 1. The above results suggested that the findings were generally robust to over-dispersion, with malnutrition showing statistical significance (p = 0.032) and ART coverage being marginally significant (p = 0.067). However, the extremely large IRR for ART coverage indicates possible model instability or data issues.

Spatial distribution and hotspot analysis of TB notifications

District-level mapping of average TB notifications over the five years showed the highest notification densities in Chipata, Katete, and Petauke districts, with annual totals ranging from 1,401 to 3,180 cases. Lusangazi and Chasefu districts recorded the lowest notification densities (<510 cases), while the remaining districts fell within the intermediate range of 410 to 1,401 cases, as shown in Figure 1. Spatial hotspot analysis using the Getis-Ord Gi* statistic identified statistically significant spatial clustering of TB notifications across Eastern Province. In 2021, TB hotspots (Z-score > 1.96) were identified in Vubwi and Chadiza districts. By 2024, significant hotspots had shifted to Sinda and Nyimba districts. Cold spots (Z-score < -1.96) were consistently identified in Chipata, Katete, and Petauke districts across the study period. Moderate clustering (Z-scores between 1.0 and 1.5) was observed in Mambwe and Chipangali districts, as shown in Figure 2.

Bayesian disease mapping results

Bayesian disease mapping using the BYM model demonstrated marked spatial heterogeneity in district-level TB risk after adjusting for population size and spatial dependence. Smoothed relative risk (RR) estimates varied across districts, with values ranging from 0.17 to 1.08. Chadiza district exhibited the highest smoothed relative risk (RR range: 0.74-1.08) as shown in Figure 3. Eight districts, Chipata, Nyimba, Katete, Lumezi, Mambwe, Vubwi, Petauke, and Kasenengwa, showed moderate relative risk estimates ranging from 0.54 to 0.74. The remaining districts, including Chama, Lundazi, Chasefu, Sinda, and Lusangazi, recorded lower relative risk estimates below 0.54. Bayesian regression estimates indicated associations between TB incidence and district-level covariates. The posterior mean coefficient for poverty was 0.4 (95% CrI: 0.1-0.7), for HIV prevalence was 0.8 (95% CrI: 0.5-1.1), and for malnutrition was 0.2 (95% CrI: -0.1-0.5) as shown in Table 3.

Discussion     

The average TB notification was 188.6 cases per district per year, cumulatively contributing to 14,145 TB notifications across the 15 districts included in this study. The observed TB notification levels are comparable to WHO estimates reported for similar predominantly rural settings. This may be attributed to limited access to TB diagnostic services, geographical barriers, and possible under-detection rather than low underlying TB transmission. The analysis revealed significant geospatial disparities in TB notifications, with Petauke, Nyimba, and Katete recording the highest TB notification counts during the study period. These districts reported the largest TB burden, which may be influenced by factors such as population size, HIV co-infection, cross-border movement, and health system access, as documented in similar settings [7]. Conversely, districts such as Chipata recorded relatively lower TB notification counts, potentially reflecting differences in healthcare access, reporting practices, or population distribution.

Malnutrition showed variability across districts; however, it was not statistically associated with TB notifications in the Poisson regression analysis, while Bayesian estimates suggested a positive but uncertain association. This mixed finding is consistent with existing evidence demonstrating a biological link between undernutrition and TB susceptibility, where malnutrition compromises immunity and increases the risk of active TB disease and poor treatment outcomes [9,10]. These results underscore the continued importance of integrating nutritional support into TB care, particularly in districts with high levels of undernutrition [11].

Spatial hotspot analysis identified localized clustering of TB notifications, with significant hotspots observed in selected districts during specific years of the study period. While Petauke, Katete, and Nyimba recorded high TB notification counts, Bayesian disease mapping further highlighted Chadiza as a district with elevated smoothed relative risk. These findings are consistent with studies from other high TB-burden settings, where spatial clustering of TB has been linked to socioeconomic and environmental factors [12,13]. The presence of both high-notification districts and areas of elevated relative risk underscores the value of combining conventional mapping with Bayesian spatial methods. Bayesian disease mapping demonstrated marked variation in district-level relative risk across Eastern Province. Chadiza exhibited the highest smoothed relative risk (0.74-1.08), while eight districts showed moderate risk levels (0.54-0.74). These spatially smoothed estimates provide a more stable representation of underlying TB risk and help identify districts that may benefit from enhanced surveillance, diagnostic capacity, and community-based interventions [11,14].

The urban-rural disparity appeared as a significant factor, with urban districts showing a 73% higher TB notification rate, equating to 88 additional cases per 100,000 population annually. HIV prevalence, though statistically significant, confirmed a clinically meaningful effect, with a 10-percentage-point increase correlating to a 45% rise in notifications, underscoring the need for integrated TB-HIV services. On the other hand, the effect of poverty was statistically significant but practically minimal (IRR = 1.00014), suggesting district-level measures inadequately capture socioeconomic distinctions. These findings highlight that urbanicity and HIV are key drivers for resource allocation, while poverty requires individual-level assessment. Overall, the practical significance of these associations varies considerably, guiding targeted interventions [15]. This analysis reaffirms the importance of distinguishing statistical from practical significance in policy decisions.

This study had several limitations. First, it relied on routine secondary data, which may be affected by reporting gaps, inconsistencies, and under-notification. Second, the ecological study design limits causal inference and may obscure individual-level variability. Finally, district-level covariates may not fully capture micro-level determinants such as household crowding, mobility patterns, and health-seeking behavior. Future studies should incorporate individual-level data, finer spatial resolution, and longitudinal designs to better understand TB transmission dynamics and validate observed spatial patterns.

Conclusion     

This study demonstrated substantial geographic variation in the tuberculosis burden across Eastern Province, Zambia, with some districts recording high notification counts and others exhibiting elevated spatially smoothed relative risk. Urban location and ART coverage were associated with higher TB notification rates at the district level. Spatial hotspot analysis and Bayesian disease mapping provided important insight into localized TB patterns and identified priority districts for intensified surveillance and targeted interventions. These findings highlight the value of integrating spatial epidemiological approaches into routine TB surveillance to support more targeted and efficient TB control efforts in hotspot areas.

Competing interests     

The authors declare no competing interests.

Authors' contributions     

Conceptualization: Claude Mumba, Aniset Kamanga. Methodology: Claude Mumba, Kaziwe Sikambale, Angel Mubanga. Data curation: Nancy Chanda-Kasese, Lindizyani Mfune, Bowa Chanda, Kabaso Mwewa. Formal analysis: Claude Mumba, Aliness Dombola, Webster Chewe, Benson Malambo Hamooya. Investigation: Claude Mumba, Chitalu Chanda, Webster Chewe. Visualization: Claude Mumba, Webster Chewe, Benson Malambo Hamooya. Writing -original draft: Claude Mumba. Writing - review and editing: Claude Mumba, Chitalu Chanda, Webster Chewe, Aniset Kamanga, Nancy Chanda-Kasese, Bowa Chanda, Kabaso Mwewa. Supervision: Aniset Kamanga, Chitalu Chanda, Benson Malambo Hamooya. All authors have read and approved the final version of this manuscript.

Acknowledgments     

The authors acknowledge the Ministry of Health, the National Tuberculosis and Leprosy Programme, and the Eastern Province Provincial Health Office for granting permission to use routine tuberculosis surveillance data. We also acknowledge the District Health Information System (DHIS2) teams for data management and reporting support.

Tables and figures     

Table 1: descriptive characteristics of districts in Eastern Province, Zambia (2020-2024)

Table 2: factors associated with TB notification rates using Poisson regression

Table 3: Bayesian spatial regression results from the Besag-York-Mollié model

Figure 1: average tuberculosis notifications in Eastern Province, Zambia (2020-2024)

Figure 2: spatial hotspot analysis using the Getis-Ord Gi* statistic for Eastern Province, Zambia (2020-2024)

Figure 3: smoothed tuberculosis relative risk estimates for Eastern Province, Zambia

Annex     

Annex 1: supplementary material (316KB)

References