Methods

Study design and analytical framework
This was a retrospective ecological analysis (2015–2023) with a hierarchical data structure examining malaria transmission in Uganda. The primary unit of analysis and statistical inference was the district-month level, comprising 13 high-burden districts observed over 108 months, yielding 1,404 district-month observations (Figure 1).

This analytical level was selected because: (1) malaria outcomes are validated and routinely reported at the district-monthly level in DHIS2; (2) individual case locations are unavailable due to privacy protections; (3) findings are directly policy-relevant for district-level planning and resource allocation; and (4) district-level aggregation provides sufficient statistical power while minimising aggregation bias.

Study area and district selection
This study was conducted in four ecologically distinct regions of Uganda: West Nile (high malaria burden, lush vegetation), Karamoja (semi-arid, sporadic transmission), Central (peri-urban landscapes), and Western (lower malaria burden, mixed land cover). Thirteen high-burden districts were selected across these regions based on the malaria caseloads (three from West Nile, three from Karamoja, four from Central, three from Western) based on: (i) the highest documented malaria caseloads during the 2018–2019 malaria indicator survey, and (ii) >80% facility reporting completeness for malaria cases among children under five years in DHIS2.

Facility and environmental data collection
Facility selection: Within the 13 study districts, eight Health Centre III facilities were selected for high-resolution geospatial environmental data collection based on: (i) highest confirmed malaria caseloads among children under five years; (ii) consistent, high-quality HMIS reporting with ≥80% parasitological confirmation rates; and (iii) geographic distribution capturing intra-district ecological diversity (approximately two facilities per district) (Table 1).

Environmental exposure measurement: Environmental variables were extracted within 5-km buffers around facility coordinates. This buffer distance was chosen to approximate facility catchment areas based on: (i) Uganda Ministry of Health guidelines for Health Centre III catchments; (ii) population mobility studies indicating 60–70% of patients reside <5 km from their primary facility; and (iii) local transport accessibility patterns  [30].

Data aggregation: Facility-level environmental measurements were aggregated to the district-month level using population-weighted averages and linked to corresponding district-month Malaria case data for all statistical analyses. Sensitivity analyses at the facility-month level validated that district-level aggregation did not substantially bias environmental effect estimates (e.g., NDVI: IRR=1.18 district-month vs. IRR=1.16 facility-month, Table 2).

Variable selection and measurement
Malaria outcome data: Malaria case data were derived from DHIS2, comprising all confirmed cases among children under five years reported by facilities within each district meeting ≥80% reporting completeness.

Environmental variables
NDVI (MODIS, 250-m resolution): Continuous in regression analyses; categorised into tertiles (low <0.3, moderate 0.3–0.6, high >0.6) for interpretability. Selected because vegetation density influences mosquito larval habitat suitability and adult vector survival.

LULC (Sentinel-2, 10-m resolution): Classification identified agriculture, shrubland, grassland, water, built-up, and bare soil. Each LULC fraction was modeled as a proportion (0–1) within facility buffers, aggregated to the district level. Selected because LULC patterns determine habitat distribution and human-environment contact.

Climate Variables (Temperature & Rainfall): Monthly minimum/maximum temperatures and rainfall were included because they directly affect parasite development and vector lifecycle completion. Variables were retained if: (i) they showed theoretical association with Malaria case, (ii) VIF <5 (collinearity threshold), and (iii) 10-year bivariate correlation with outcome |r| >0.10.

Justification for Multi-Scale Approach: Fine-scale facility-level exposure measurement (5-km buffers) was followed by aggregation to the district-month analytical unit, leveraging high-spatial-resolution environmental data while aligning with outcome data structure and policy relevance. Facility-level GIS spatial analysis verified that vegetation-malaria associations observed at the district-month level reflect mechanisms operating at fine spatial scales.

Data management and missing data handling
Malaria Surveillance Data: Missing values occurred in 4.2% of district-month observations, primarily due to delays in routine reporting. These were handled using multiple imputation by chained equations (MICE) under the missing-at-random (MAR) assumption. Ten imputed datasets were generated using auxiliary variables (district population size, health facility density, and seasonal trends). Regression models were run across imputed datasets with results pooled using Rubin’s rules. Complete-case analysis was performed as a sensitivity analysis to assess robustness.

Environmental Data: Missing values were minimal (<1%) and addressed through linear temporal interpolation between adjacent months. Spatial interpolation using inverse distance weighting (IDW) from neighboring districts was applied when necessary. All imputed values were cross-checked against historical seasonal norms.

Analytical unit and data aggregation
This study employed the district-month as the primary unit of analysis and statistical inference, comprising 13 high-burden districts observed over 108 months (2015–2023), yielding 1,404 district-month observations. Each observation represented Malaria case data and environmental exposures for one district in one calendar month.

Justification for district-month level: The district-month level was selected for four key reasons: (1) malaria outcomes in DHIS2 are validated and routinely reported at this temporal and spatial resolution; (2) individual case locations are unavailable due to privacy protections; (3) findings are directly actionable for district-level planning and resource allocation; and (4) this aggregation provides sufficient statistical power for robust inference while minimizing aggregation bias.

Data aggregation process: Malaria case data were derived from all health facilities within each district that met >80% reporting completeness criteria. Eight Sentinel Health Centre III facilities (approximately two per district) were selected based on the highest malaria case volumes to guide environmental data collection. Environmental variables (NDVI, LULC, rainfall, and temperature) were extracted within 5-km buffers around these facilities to characterize district-level exposures. Facility-level measurements were subsequently aggregated to the district-month level using population-weighted averages and linked to corresponding Malaria case data for all statistical analyses.

Multi-scale approach: This methodological design leveraged high-resolution satellite-derived environmental data (NDVI at 250-m resolution, LULC at 10-m resolution) at the facility scale while maintaining the district-month as the analytical unit. Sensitivity analyses at the facility-month level (864 observations) confirmed that district-level aggregation did not substantially bias environmental effect estimates (e.g., NDVI: IRR=1.18 at district-month vs. IRR=1.16 at facility-month), validating this aggregation strategy. All primary results and inferences are reported at the district-month level; no facility-level regression analyses were performed.

Data sources
Malaria case data
Malaria case counts for children under five years were obtained from Uganda’s District Health Information System version 2 (DHIS2) for 2015–2023. Data were cleaned to remove duplicates and missing values, then aggregated into monthly district-level totals to construct continuous time series enabling analysis of seasonal, inter-annual, and lagged relationships with environmental covariates.

Environmental data: Specifications and harmonisation
Vegetation (NDVI): Monthly mean NDVI was derived from the MODIS MOD13Q1 product (250-m resolution, 16-day intervals). Within-month observations were averaged to align with monthly malaria data, reducing short-term variability while preserving seasonal patterns. Missing values (<1%) due to cloud cover were addressed through linear interpolation.

Land Use and Land Cover (LULC): Annual LULC maps were generated from Sentinel-2 multispectral imagery (10-m resolution) using supervised classification in Google Earth Engine, classified into five categories: trees, shrubs, agriculture, wetlands, and built-up areas. District-level proportions were calculated by overlaying administrative boundaries on annual maps. Due to the slow temporal dynamics of land cover, annual LULC proportions were uniformly assigned to all 12 months within each calendar year, a standard approach in malaria-environment research. This temporal approximation was balanced by incorporating monthly NDVI as a dynamic vegetation indicator and applying monthly fixed effects to account for seasonal patterns Table 3.

Climate variables: Monthly rainfall was retrieved from CHIRPS v2 (0.05° resolution, ~5 km). Minimum and maximum 2-meter air temperatures were obtained from ERA5 reanalysis (0.1° resolution, ~9–11 km). All values were extracted and averaged within 5-km buffers around selected health facilities. Table 3 shows the summary of variables and their measurements

Dataset integration
Each monthly district-level observation contained: (1) one malaria case count, (2) one NDVI value, and (3) eight unchanging LULC proportions (summing to 100%) corresponding to that calendar year. This structure enabled panel time-series modelling while leveraging high-resolution satellite data and aligning with malaria reporting frequencies.

Temporal and statistical analysis
Temporal decomposition
Seasonal-Trend decomposition using LOESS (STL) was applied to Malaria case, NDVI, and LULC time series to separate trend, seasonal, and residual components while accommodating non-linear trends and variable seasonality. A 12-month seasonal window with robust smoothing was implemented in R to isolate meaningful long-term patterns from short-term fluctuations.

Environmental exposure measurement and linkage
Environmental variables were extracted at the facility level within 5-km catchment buffers, which align with Uganda Ministry of Health guidelines for Health Centre III service areas and epidemiological evidence of patient travel patterns. These facility-level measurements were then aggregated to the district-month level, the primary unit of analysis, to match malaria outcome reporting structures and achieve sufficient statistical power.

NDVI (MODIS, 250-m resolution) was averaged monthly to align with monthly Malaria case data. For subgroup analyses examining dose-response associations, NDVI was categorized into tertiles (low <0.3, moderate 0.3–0.6, high >0.6) to enhance the interpretability of vegetation effects.

LULC (Sentinel-2, 10-m resolution) was mapped annually, producing eight land-cover proportions per district per year (trees, shrubs, agriculture, and wetlands, built-up areas, grassland, water, and bare soil). Due to the slow rate of land-cover change relative to monthly malaria fluctuations, these annual proportions were assigned uniformly to all 12 months within each calendar year, a standard approximation in malaria-climate research. Complementary monthly NDVI data and monthly fixed effects in regression models mitigated potential temporal misalignment in highly dynamic landscapes.

Climate variables were extracted monthly: rainfall from CHIRPS (0.05° resolution, ~5 km) and minimum/maximum temperatures from ERA5 reanalysis (0.1° resolution, ~9–11 km). These were averaged within facility buffers to match the spatial scale of Malaria case and vegetation data.

This multi-scale approach, fine-resolution facility-level extraction followed by district-month aggregation, enabled the use of high-resolution satellite data while maintaining consistency with routinely reported health outcomes and policy-relevant geographic units.

Statistical analysis
Correlation analysis: Pearson’s correlation coefficients quantified linear relationships between Malaria case (outcome) and environmental/climatic predictors. Vegetation cover was operationalized as the combined proportion of trees, shrubs, and agricultural land, providing a broad measure of vegetated habitat. Prior to analysis, variables were examined for distributional properties and normality; transformations were applied where necessary to meet Pearson correlation assumptions. Statistical significance was evaluated at the 95% confidence level (α = 0.05) with notation: *** (p < 0.001), ** (p < 0.01), * (p < 0.05).

Regression modelling: Fixed-effects negative binomial regression models were used to examine associations between aggregated environmental exposures and malaria cases among children under five years. This modeling approach accounts for the count nature of malaria data and the wide range of baseline incidence across districts. Standard errors were clustered at the district level to account for within-district temporal correlation. Incidence Rate Ratios (IRRs) with 95% confidence intervals are reported, each interpreted as the relative percentage change in Malaria cases associated with a one-unit increase in the corresponding environmental predictor.

Diagnostic testing and robustness: Temporal autocorrelation was assessed using the Wooldridge test, which confirmed significant serial correlation in the data. To address this, two-way cluster-robust standard errors were estimated at the district level to account for both within-district correlation and temporal dependence. Spatial dependence was evaluated using Moran’s I test on model residuals. Sensitivity analyses using Newey-West heteroskedasticity-autocorrelation consistent standard errors confirmed the robustness of effect estimates. Post-estimation diagnostics showed minimal residual autocorrelation, with finite-sample bias corrections applied to variance-covariance matrices to ensure valid inference.

Lag structure
Environmental covariates were evaluated both at the same time (lag 0) and with 1-month (lag 1) and 2-month (lag 2) delays relative to reported malaria cases. This lagged approach captures delayed effects on transmission caused by biological and operational factors. Specifically, rainfall and suitable vegetation increase breeding habitats for Anopheles mosquitoes, leading to rises in adult mosquito populations roughly 1–2 weeks after habitat improvements. After mosquito development and infected blood meal intake, the Plasmodium parasite experiences an extrinsic incubation period of about 10–14 days (temperature-dependent), further delaying infectious transmission. Additionally, delays in surveillance reporting within routine health systems can shift observed malaria cases relative to actual environmental conditions. Together, these biological and operational factors support including 1–2 month lag structures to assess both immediate and delayed links between environment and malaria.

Spatial visualization and software
Geographic heterogeneity in Malaria case and environmental covariates was visualized at the health facility level using heat maps generated in ArcGIS Pro to identify spatial patterns. However, primary analyses focused on correlation and regression modeling to rigorously quantify associations between environmental factors and malaria outcomes. All statistical analyses were conducted in R (version 4.2.1) to ensure transparency and reproducibility of results

Model selection and justification
Given the count nature of monthly malaria cases and evidence of overdispersion (variance exceeding the mean), we employed a negative binomial panel regression model. The model specification was:

where:

• $\mu_{it}$ is the expected malaria case count for district $i$ at month $t$
• $\log(\text{Population}_{it})$ is the offset term, the natural logarithm of the under-five population included to model incidence rates rather than raw counts (coefficient constrained to 1)
• $u_i$ are district-specific random effects, assumed $u_i \sim N(0, \sigma_u^2)$
• $\varepsilon_{it}$ is the error term.

The negative binomial distribution accommodates overdispersion through the dispersion parameter $\alpha$, where $\text{Var}(Y_{it}) = \mu_{it} + \alpha\mu_{it}^2$. A statistically significant $\alpha > 0$ confirms overdispersion and justifies the negative binomial over Poisson specification.

To account for unobserved heterogeneity across districts, both fixed-effects and random-effects negative binomial models were estimated. A Hausman test was conducted to choose between the fixed effects and random effects models. The test statistic was χ² (4) = 25.7118 (p = 0.0000), leading to the rejection of the null hypothesis. Therefore, the fixed effects model was selected for the analysis. Consequently, the fixed-effects negative binomial model was adopted as the primary analytical model.

Model diagnostics further supported this choice in Table 5. The model showed a good overall fit. Model performance metrics for the fixed-effects specification were: Log Likelihood = −5852.0, Akaike Information Criterion (AIC) = 11740.0, and Bayesian Information Criterion (BIC) = 11832.3. The dispersion parameter (alpha = 0.1434) confirmed the appropriateness of the negative binomial model over Poisson regression. Residual diagnostics showed a residual mean close to zero (0.0017) and acceptable distribution characteristics (skewness = 1.3374, kurtosis = 3.4370). Nevertheless, random-effects estimates are also presented as a robustness check because they allow estimation of between-district variation and facilitate comparison of effect sizes across model specifications. Comparing results across both approaches provides an assessment of the consistency and robustness of the observed relationships between climate variability and Malaria cases.

Incorporation of temporal structure and handling of spatial dependence
To adequately account for the strong temporal structure characteristic of malaria transmission, monthly dummy variables were included to capture seasonality driven by rainfall and temperature patterns. As a robustness check, model seasonal cycles were also tested. Potential lag effects of environmental covariates were examined by including 1-month and 2-month lagged values, selected based on exploratory cross-correlation analysis and model information criteria. A dynamic specification with the lagged dependent variable was additionally explored where appropriate.

Spatial dependence and within-district temporal autocorrelation were addressed in multiple ways. First, district-level fixed effects absorbed time-invariant spatial heterogeneity. Second, standard errors were clustered at the district level to allow for arbitrary within-district correlation over time (including serial correlation). The presence of serial correlation was formally tested using Wooldridge’s test for autocorrelation in panel data.

Confounding and collinearity
Confounding was addressed using multiple complementary strategies. To control all unobserved time-invariant district-level confounders (e.g., topography, socio-economic conditions, health infrastructure, and intervention coverage), district fixed effects were included in the negative binomial regression models. Time-varying confounders were adjusted for by incorporating monthly mean rainfall, temperature, and relative humidity, each included with 1- and 2-month lags. Seasonality was further controlled using monthly dummy variables. Normalized Difference Vegetation Index (NDVI) was included as a dynamic, time-varying covariate to account for changes in vegetation cover.

Assessment of collinearity was performed prior to model fitting. Pairwise Pearson correlations were examined among all continuous environmental predictors (rainfall, temperature, humidity, and NDVI). In addition, Variance Inflation Factors (VIF) were calculated for each covariate in the final multivariable models. All VIF values were below 5 (mean VIF < 2.5), indicating that multicollinearity was not a concern. Where moderate correlations were observed (e.g., between temperature and humidity), models were tested with and without one of the correlated variables to confirm the robustness of the primary associations.

Ethical approval
Ethical approval for this study was obtained from the Mulago Hospital Research and Ethics Committee (MHREC-2910). Administrative clearance and authorization to access and use the data were granted by the Ministry of Health, Uganda. Permission was specifically granted for the use of secondary surveillance data extracted from Uganda’s national Health Management Information System (HMIS/DHIS2).

Given that the study relied exclusively on routinely collected, aggregated secondary surveillance data with no direct contact with human participants, an exemption from informed consent (waiver of consent) was granted by the ethics committee. All data were fully anonymized at the source before access, and no individual-level identifiers were obtained or analyzed at any stage of the study.

Data extraction from the HMIS was conducted using a standardized data abstraction tool developed for this study to ensure consistency in the retrieval of malaria case data across facilities and time periods. The tool facilitated structured extraction of aggregated monthly malaria reports and relevant metadata for analysis.

All procedures were conducted in accordance with national research ethics guidelines and the principles of the Declaration of Helsinki. Confidentiality and data security were strictly maintained throughout the study, and all analyses were performed on de-identified aggregated datasets.