The influence of malaria control interventions and climate variability on changes in the geographical distribution of parasite prevalence in Kenya between 2015 and 2020

Abstract

Background

The burden of malaria in Kenya was showing a declining trend, but appears to have reached a plateau in recent years. This study estimated changes in the geographical distribution of malaria parasite risk in the country between the years 2015 and 2020, and quantified the contribution of malaria control interventions and climatic/ environmental factors to these changes.

Methods

Bayesian geostatistical models were used to analyse the Kenyan 2015 and 2020 Malaria Indicator Survey (MIS) data. Bivariate models were fitted to identify the most important control intervention indicators and climatic/environmental predictors of parasitaemia risk by age groups (6–59 months and 5–14 years). Parasitaemia risk and the number of infected children were predicted over a 1 × 1 km² grid. The probability of the decline in parasitaemia risk in 2020 compared to 2015 was also evaluated over the gridded surface and factors associated with changes in parasitaemia risk between the two surveys were evaluated.

Results

There was a significant decline in the coverage of most malaria indicators related to Insecticide Treated Nets (ITN) and Artemisinin Combination Therapies (ACT) interventions. Overall, there was a 31% and 26% reduction in malaria prevalence among children aged < 5 and 5–14 years, respectively. Among younger children, the highest reduction (50%) and increase (41%) were in the low-risk and semi-arid epi zones, respectively; while among older children there was increased risk in both the low-risk (83%) and semi-arid (100%) epi zones. Increase in nightlights and the proportion of individuals using ITNs in 2020 were associated with reduced parasitaemia risk.

Conclusion

Increased nightlights and ITN use could have led to the reduction in parasitaemia risk. However, the reduction is heterogeneous and there was increased risk in northern Kenya. Taken together, these results suggest that constant surveillance and re-evaluation of parasite and vector control measures in areas with increased transmission is necessary. The methods used in this analysis can be employed in other settings.

Introduction

Malaria causes substantial disease and deaths, especially in sub-Saharan Africa (SSA). It was estimated that in 2022, there were approximately 249 million cases and 608,000 deaths worldwide [1]. The World Health Organization (WHO) Afro-region contributed 94% of cases and 96% of deaths due to malaria. Efforts have been made to accelerate the reduction of malaria burden. In 2015, the World Health Assembly (WHA) adopted the Global Technical Strategy (GTS) for Malaria (2016–2030) to reduce malaria incidence and deaths by at least 90% by 2030 [2]. However, in 2017, the GTS’s aim was off-course, leading to the launch of the High Burden to High Impact (HBHI) approach in 2018 by the WHO and the Roll Back Malaria partnership [3]. The HBHI approach supports 11 countries that account for approximately 70% of the global burden to reduce their malaria burden and achieve the GTS 2025 milestones [3].

Similar strategies have been adopted and implemented in Africa. For instance, the Africa Union (AU) initiated the “Zero Malaria Starts with Me” campaign in 2021 with the goal of eliminating malaria by 2030 [4]. Additionally, individual countries have developed strategies to fight malaria. For example, Kenya introduced the “Kenya Malaria Strategy (KMS) 2019–2023” which aimed at reducing malaria incidence and deaths by 75% of the 2016 levels by 2023 [5]. The success of these strategies depends upon several factors including; optimal use of effective vector control tools, variability in climatic and environmental factors, community awareness and involvement, strengthening the capacity of National Malaria Control Programs (NMCP) and encouraging partnerships [2, 5,6,7,8]. Monitoring the spatio-temporal distribution of parasitaemia risk, and estimating the effectiveness of interventions are essential in guiding implementation of targeted control measures and resource allocation.

Several studies have used geostatistical methods to map country-specific malaria risk in SSA [9,10,11,12,13,14,15,16,17,18] with the aim of informing NMCPs. In Kenya, earlier studies mapped malaria risk using empirical methods based on parasitological, demographic, geographical and climatic data [19, 20]. However, more recent studies have employed model-based geostatistical mapping [10, 14, 17]. These maps provide important information to the Kenyan NMCP about malaria zones and spatio-temporal changes in malaria risk over time, offering valuable insights for malaria control.

Mategula and Gichuki evaluated the heterogeneity and spatial drivers of malaria transmission at a 5 × 5 km² spatial resolution, using the Kenya 2020 malaria indicator survey (MIS) data [21]. In their analysis, they aggregated the data across all ages for children aged 0.5 to 14 years. However, the risk and severity of malaria differs by age group with children aged < 5 years more likely to have severe outcomes [22], while school aged children (5–14 years) act as reservoirs of malaria parasites [23] without showing symptoms; therefore the probability of being tested may have been different. Additionally, most studies assess malaria risk in children < 5 years old, where control interventions are targeted, however, there is evidence to show that the risk of malaria shifted to older children after scale up of interventions [24]. Therefore, it is necessary to evaluate changes in malaria risk in older children as well. Furthermore, it is important to stratify the analyses by age groups, as we have previously found that aggregating ages may mask the effects of climatic predictors [25]. In addition to risk maps, it is important to inform the NMCP about the probabilities of risk reduction/ increase, and the risk by different age groups at fine spatial scales. Additionally, the joint effects of changes in control intervention coverage and climatic/environmental factors on changes in parasitaemia risk remains underexplored.

Most studies, including previous Kenyan studies that use data from different sources such as MIS, use binomial models to map parasitaemia risk [10, 14, 17, 21]. However, malaria parasitaemia prevalence in Kenya is below 10%, and most of the sampled MIS locations have zero prevalence [26, 27]. As Kenya and other countries with lower prevalence move towards malaria elimination, the commonly used binomial distribution may under-estimate the probability of zero prevalence. Although not widely adopted for prevalence data, zero-inflated binomial (ZIB) models can be used to account for the excess zeros [12].

The current study uses Bayesian geostatistical binomial and ZIB models to assess the net effect of control interventions and climatic/ environmental factors on malaria parasitaemia risk by age group, and to predict parasitaemia risk at a spatial resolution of 1 km² in Kenya during 2015 and 2020. In addition, this study estimates the population adjusted parasitaemia prevalence, maps fine scale probability of risk reduction between the two surveys, estimates number of infected children at every pixel during the two surveys and the effects of changes in control interventions and climatic/environmental factors between 2015 and 2020 on the geographical distribution of the changes in parasitaemia risk.

Materials and methods

Study setting and population

Kenya, located in East Africa, is home to 52 million people [28], has a total landmass of 580,367 square kilometres. It shares its borders with five countries: Ethiopia, Somalia, Tanzania, Uganda, and South Sudan. Kenya’s geography is diverse, featuring a long coastline along the Indian Ocean, coastal lowlands, savannahs, arid and semi-arid regions in the north, vast forests, highlands, and mountainous areas. Additionally, Kenya is divided into 47 administrative level 1 units referred to as counties with significant variability in climate, ranging from the tropical coast to the temperate inland and arid northern regions. The coolest months are from June through August, while the warmest months are from December through March. Temperature and rainfall vary with altitude, with the coastal areas, eastern plateaus, highlands, and northern plains recording temperatures between 18 and 34 °C, 14–29 °C, 10–26 °C, and 33–38 °C, respectively. Rainfall exhibits a bi-modal pattern, with long rains occurring from March to May and short rains from October to December. Total precipitation ranges from under 250 mm in the north and north eastern regions to 1,270 mm in the coastal regions.

Malaria accounts for 13–15% of outpatient consultations, with approximately 70% of the population at risk [27]. Kenya is divided into five malaria epidemiological zones, including, the highland epidemic-prone areas, lake endemic areas, coast endemic areas, semi-arid seasonal epidemic areas, and low-risk epidemic areas [5]. Plasmodium falciparum parasite is responsible for 96–99% of infections [27]. A map depicting the county outlines and the malaria epidemiological zones is provided in additional file 1, Figure A1.

The Ministry of Health distributed Insecticide-Treated Nets (ITNs) to high malaria transmission zones beginning in 2006, targeting children aged under 5 years. Subsequent mass distribution of ITNs in 2011/2012, 2014/2015 and 2017/2018 targeted all individuals living in the malaria-prone zones with the objective of achieving universal coverage (one net for two people). Targeted indoor residual spraying (IRS) for vector control was employed in 12 counties from 2006 then expanded to 16 epidemic prone and 4 endemic counties during 2010/2011 before being suspended after 2012 to address insecticide resistence issues [14]. However, in 2017, IRS activities resumed in two counties, Migori and Homa Bay, located in the Lake endemic zone [27]. Artemisinin-based Combination Therapy (ACT) and rapid diagnostic tests (RDT) were introduced in 2006 and 2012 respectively [14].

Data sources

MIS data

The MIS, designed to follow the Roll Back Malaria’s (RBM) monitoring and evaluation working group guidelines and national malaria strategies, provides representative estimates of key malaria indicators. The 2015 MIS was conducted in 6,481 households across 246 clusters between July and August, 2015; while the 2020 MIS took place in 7,952 households within 301 clusters between November and December, 2020 [26, 27]. All women aged 15–49 years in these households were eligible for interviews, and finger or heal prick blood samples were collected from children aged 6 months to 14 years for malaria testing. Children aged 6–59 months are referred to as children aged < 5 years in the text. The blood samples were tested on-site using the SD Bioline (Abbott Diagnostics, Korea) and CareStart Pf (Access Bio inc, NJ, USA) RDTs in 2015 and 2020, respectively, with results provided to the participants. Children who tested positive by RDT were either treated or referred to a heath facility for immediate treatment. Additionally, thick and thin blood films were taken to malaria reference laboratories for microscopy testing. For this analysis, we used data from 245 clusters in 2015 and 298 clusters in 2020, considering malaria positivity by microscopy. The cluster locations are displayed in Figs. 1 and 2.

Fig. 1

Observed malaria prevalence in children aged < 5 years at survey locations in 2015 (A) and 2020 (B); predicted parasitaemia prevalence in 2015 (C) and 2020 (D); coefficient of variation of the predicted prevalence in 2015 (E) and 2020 (F); estimated number of infected children per km² in 2015 (G) and 2020 (H)

Fig. 2

Observed malaria prevalence in children aged 5–14 years at survey locations in 2015 (A) and 2020 (B); predicted parasitaemia prevalence in 2015 (C) and 2020 (D); coefficient of variation of the predicted prevalence in 2015 (E) and 2020 (F); estimated number of infected children per km² in 2015 (G) and 2020 (H)

Enviromental/ climatic data

Downscaled monthly day Land Surface Temperature (LSTD), night Land Surface Temperature (LSTN) and Normalized Difference Vegetation Index (NDVI) data were extracted from the Moderate Resolution Imaging Spectroradiometer (MODIS) [29]. We obtained downscaled rainfall data from the Climate Hazards Group Infrared Precipitation with Station data (CHIRPS) [30]. Nightlights – which is a proxy for urbanization and socio-economic status (SES) [31] – was obtained from National Oceanic and Atmospheric Administration - Visible Infrared Imaging Radiometer Suite (NOAA-VIIRS) [32] while land cover data was obtained from Copernicus [33]. Distance to permanent water bodies was estimated as the shortest Euclidean distance from the cluster centroids to the nearest major water body using the 100 × 100 m land cover data from Copernicus. Altitude for each cluster was estimated from the Shuttle Radar Topographic Mission (SRTM) [34] using the digital elevation model. The climatic indicators of temperature, vegetation, and rainfall were summarised by calculating their average values over the six-month period preceeding the end of the MIS. These average values were extracted at the survey cluster locations. Detailed descriptions of these data sources and their spatial and temporal resolutions are displayed in additional file 1, Table A1.

Interventions

Data on ITN ownership and use, ACT use, and health seeking were extracted from the MIS. Based on these data, indicators were calculated using the definitions provided by RBM [35]. These indicators included: (i) the proportion of households with at least one ITN; (ii) the proportion of households with one ITN for every two people; (iii) the proportion of the population with access to an ITN within their household; (iv) the proportion of the population that slept under an ITN the previous night; (v) the proportion of children under 5 years old who slept under an ITN the previous night (vi) the proportion of existing ITNs used the previous night; (vii) the proportion of children with fever who sought care; and (viii) the proportion of fever episodes treated with ACTs. Notably, we excluded IRS use from this analysis as this variable was not collected during the 2020 MIS.

Population data

Gridded high-resolution population data for 2015 and 2020 were obtained from Worldpop at a 100 × 100 m² spatial resolution [36]. The United Nations Population Division data on population structure [37] were applied to calculate the number of children aged < 5 and 5–14 years old for each year. Specifically, the Worldpop data were multiplied by the proportion of children in each age group. For children aged < 5 years, the proportions were 14.5% and 13.2% in 2015 and 2020, respectively, while for those aged 5–14 years, the proportions were 27.0% and 25.7% in 2015 and 2020, respectively.

Statistical analyses

Data management and statistical analyses were conducted in R version 4.1.3 [38]. Geostatistical model and variable selection, model fitting and kriging were carried out using Integrated Nested Laplace Approximation (INLA) [39]. Climatic/environmental variables were categorized into groups based on their tertiles, which divides the data into three equal parts. Both continuous and categorical forms of the variables were assessed for their association with parasitaemia risk. The selection criterion was based on the Deviance Information Criterion (DIC). Priority was given to the form with the lowest DIC, except when the difference in DIC was smaller than 5, in which case the continuous form was prioritized. Variables were standardized (scaled using their mean and standard deviation) to allow for a relative comparison of covariate effects.

Bayesian geostatistical ZIB models were fitted to the MIS data separately at each survey year and age group through stochastic partial differential equations (SPDE) to: (i) estimate the effects of control interventions and climatic/environmental factors on parasitaemia risk; (ii) obtain spatially explicit predictions of parasitaemia prevalence and of the probability of a decline between the two survey years; (iii) estimate the population-adjusted prevalence, number of children infected with malaria parasites and their relative reduction by epidemiological zone and county; and (iv) assess the effects of the changes in coverage of control intervention coverage and climatic factors between 2015 and 2020 on the geographical distribution of the changes in parasite risk. Bivariate geostatistical models were fitted to assess the association between individual predictors and parasitaemia risk and to determine the best distributional assumption between binomial and ZIB based on the DIC. Rainfall and nightlights were included in the multivariable models a priori guided by previous work in Africa [10]. Subsequently, other covariates including socio-demographic, climatic and intervention factors were added and the model’s fit and predictive ability evaluated.

Bayesian model formulation

Let $\:{Y}_{i}$ be the number of children who tested positive, and $\:{N}_{i},$ be the total number of children tested in cluster $\:{\text{s}}_{\text{i}}$. Typically, $\:{Y}_{i}$ is assumed to come from a binomial distribution ${Y_i} \sim Bin({N_i},{p_i})$ where $\:{p}_{i}$ is the probability of parasitaemia at cluster $\:{\text{s}}_{\text{i}}.$ However, in cases where excess zeros are observed, the binomial distribution may not be able to estimate the zero prevalence probability or correctly identify covariates associated with the outcome. We fitted a ZIB model ${Y_i} \sim ZIB\left( {{N_i},{p_i},{\theta _i}} \right)$ which assumes two sources of zeros: $\:{\theta\:}_{i}$% (mixing probabilities) of structural zeros and $\:\left({1-\theta\:}_{i}\right)\%$ arising from the binomial distribution. The model is defined as

$$\:{Y}_{i}|{p}_{i},{\theta\:}_{i}\sim\left\{\begin{array}{c}0\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:with\:probability\:{\theta\:}_{i}\\\:Bin\left({N}_{i}\:,{p}_{i}\:\right)\:\:\:\:\:\:\:with\:probability\:1-{\theta\:}_{i}\end{array}\right.$$

where the probability $\:{\theta\:}_{i}$ is modelled using covariates and is defined as $\:{logit(\theta\:}_{i})=\sum\:_{k=1}^{m}{\alpha\:}^{T}\varvec{X}\left({s}_{i}\right)$ where $\:\varvec{X}\left({s}_{i}\right)$ are covariates at location ${s}_{i}$ and $\:\alpha\:$ = ($\:{\alpha\:}_{1},\:{\alpha\:}_{2},\dots\:,{\alpha\:}_{k}$)^T are the regression coefficients. The relationship between $\:{p}_{i}$ and the predictors was modeled via the equation $\:logit\left({p}_{i}\right)={\beta\:}^{T}\varvec{X}\left({s}_{i}\right)+\omega\:\left({s}_{i}\right)$ where $\:\varvec{X}\left({s}_{i}\right)$ is the set of predictors at location$\:\:{s}_{i}$, $\:\varvec{\beta\:}$ = ($\:{\beta\:}_{1},\:{\beta\:}_{2},\dots\:,{\beta\:}_{k}$)^T is the vector of regression coefficients and $\:\omega\:\:$are the spatial random effects included in the model via the Matern covariance function. The ZIB models were fitted using INLA’s 0binomial function.

A regular grid at 1 × 1 km² resolution was created, and parasitaemia risk was predicted at the centroids of the grid cells (unsampled locations) through Bayesian kriging based on the model with the best predictive ability (the least root mean square error (RMSE), the least mean absolute error (MAE) and the highest correlation coefficient at each survey year and age group). For each age group, the probability of reduction in parasitaemia prevalence in 2020 compared to 2015 was estimated. In particular, 100 samples were drawn from the posterior predictive distribution of the parasitaemia prevalence at each grid cell. The difference in prevalence between these samples was calculated, with a value of one assigned if there was a decline, and zero otherwise. These sample-based differences were then averaged to obtain the probability of parasitaemia reduction. The posterior predictive distribution of the number of infected children was determined by multiplying the posterior samples at each grid-cell with the respective, age-specific population data. The overall, epidemiological zone-specific, and county-specific numbers of infected children were calculated by summing the pixel level estimates within the country, zones and counties, respectively. The population adjusted prevalence was obtained by dividing the total number of infected children by the age-specific population. Quantiles of the sample-based posterior predictive distribution of the population-adjusted prevalence were used to calculate 95% Bayesian credible intervals (BCIs) for the above quantities.

A separate analysis was conducted to assess the effects of the changes in coverage of control intervention indicators and climatic factors between 2015 and 2020 on the geographical distribution of the changes in parasitaemia risk. In particular, a Bayesian geostatistical ZIB model was fitted to the parasitaemia survey data of 2020. The parasite prevalence of 2020 at a given location, transformed on the logit scale, was modelled as a function of the logit transformed parasitaemia prevalence of 2015 (with the regression coefficient constrained to be equal to 1), the effect of control intervention coverage and climatic/environmental factors in 2020 adjusted for the difference in control interventions coverage and the difference of climatic/environmental factors between the two survey years. Geographical misalignment of the cluster locations between the two surveys was addressed by extracting the predicted parasitaemia risk of the 2015 survey and 2015 climatic covariates at the 2020 survey locations. The difference in ITN coverage between the two surveys was estimated at the county level. To increase the number of survey locations, we included the 2014 demographic and health survey (DHS) clusters to the 2015 MIS clusters to estimate the ITN coverage during 2015 [26, 40].

Odds ratios (OR) and adjusted ORs (aOR), quantifying the effect of the predictors were summarized by their posterior median and the corresponding 95% BCIs. The effects were considered statistically important if their 95% BCI did not include one. Throughout this paper, we adhere to terminology consistent with Bayesian inference. Specifically, odds ratios that exclude one in their 95% BCIs are referred to as statistically important, rather than “statistically significant” commonly used in frequentist inference.

Ethical consideration

The protocols for the Kenya Malaria Indicator Surveys were reviewed and approved as previously described [26, 27]. We use aggregated secondary data for this analysis and do not attempt to identify the individuals involved in the surveys. Further, permission to use these data for our study was sought and approval obtained from the DHS program.

Results

Descriptive analysis

This analysis included 10,037 children in 2015, with 3,429 (34.2%) aged < 5 years. In 2020, the analysis involved 11,491 children, with 3,725 (32.4%) falling in the < 5 years age group. Overall, there was a significant decline in the prevalence of parasitaemia from 8.2% in 2015 to 5.6% in 2020 (p-value < 0.01). Over the two surveys, there was a lower prevalence among children aged < 5 years compared to children aged 5–14 years; with rates of 5.0% vs. 9.9% in 2015 and 3.0% vs. 6.9% in 2020. The lake endemic epidemiological zone had the highest prevalence (26.6%) in 2015, but in 2020 there was a reduction in parasitaemia in the lake, coast-endemic, and highlands zones. In contrast, there was an increase observed in the semi-arid (from 0.5 to 1.8%) and the low-risk (from 0.3 to 0.4%) epidemiological zones. These results are summarized in Table 1 and Table A2 in additional file 1.

Table 1 Description of the 2015 and 2020 malaria indicator survey data and climatic/environmental factors

There was a decline in ITN ownership and use between the two surveys. For the ITN ownership indicators, the proportion of households with at least one ITN, the proportion of households with at least one ITN for every two people, and the proportion of individuals with access to an ITN in their household declined by 22%, 28% and 15%, respectively (Table 1). For the ITN use indicators, the proportion of individuals that reported having slept under an ITN the previous night and the proportion of children under 5 years old who slept under an ITN declined by 27% and 25%, respectively, whereas the proportion of existing ITNs used the previous night increased by 7%. Furthermore, the proportion of children with fever for whom care was sought, and the proportion of fever episodes treated with ACTs declined by 11% and 26%, respectively.

Table 1 summarizes the changes in climatic conditions. The mean rainfall, crop cover, and nightlights all increased in 2020, while there was a decrease in LST day and the distance from cluster centroids to the permanent water bodies. Fine spatial scale changes in rainfall, LSTD, nightlights and ITN use are presented in Figure A2 in additional file 1, where we observed increased rainfall in the western part of Kenya, reduced temperature in Turkana county and parts of central Kenya, and increased nightlights in all the counties, while changes in ITN use showed no clear clustering.

Model and variable selection

The geostatistical ZIB models consistently outperformed their binomial counterparts, evidenced by lower DIC values in all bivariate models. For children under 5 years old, important predictors of parasitaemia risk included rainfall, LST, nightlights, altitude, distance to permanent water bodies, and the proportion of existing ITNs used during the previous night, or the proportion of individuals that used ITN the previous night were important predictors of parasitaemia risk. Similarly, these variables, along with ITN ownership indicators, were important predictors of parasitaemia risk among older children. Detailed results can be found in Table A3 and Table A4 (Additional file 1).

For prediction of parasitaemia risk in 2015, 57 and 120 model combinations were fitted for children aged < 5 and 5–14 years, respectively. For the 2020 survey 57 and 26 model combinations were fitted for children aged < 5 and 5–14 years, respectively. The top five models for each group had comparable predictive ability based on RMSE, MAE and correlation coefficients, Table A5 (Additional file 1). The champion model in each group was used in plotting the final maps.

Effects of climatic/environmental factors and interventions on parasitaemia prevalence

Posterior estimates of the effects of control interventions and climatic factors on parasitaemia risk in 2015 and 2020 for children aged < 5 years and 5–14 years are shown in Tables 2 and 3, respectively. In 2015, increase in nightlights and altitude was associated with a reduced risk among younger children, while an increase in rainfall and the proportion of households with at least one ITN were associated with an increased risk. Nightlights was also associated with reduced risk in older children. In 2020, increase in the proportion of individuals sleeping under ITNs was associated with a reduced risk in both age groups, while an increase in the use of ACT had a positive effect.

Table 2 Effects of climatic/ environmental factors and interventions on parasitaemia prevalence in children aged 6–59 months

Table 3 Effects of climatic/ environmental factors and interventions on parasitaemia prevalence in children aged 5–14 years

Spatio-temporal trends of parasitaemia risk during 2015 to 2020

The geographical distribution of the parasitaemia risk, predicted over the 1 × 1 km² grids among children aged < 5 years and children aged 5–14 years is presented in Figs. 1 and 2, respectively. The risk among children aged < 5 years was high in the western and coastal parts of the country in both surveys. These regions had the lowest coefficient of variation, while Nairobi had the highest. Overall, the population adjusted risk reduction among younger children in 2020 relative to 2015 was 31.3%. Notably, parasitaemia risk expanded to new areas north of the country by 2020 contributing to the 41.2% increase in parasitaemia risk observed in the semi-arid seasonal zones. When compared by counties, there was a 74% reduction in parasitaemia risk in Homa Bay County, and an increase in Turkana County (3–4 times) between the two surveys.

A similar pattern was observed in the geographical distribution for children aged 5–14 years, with a risk reduction of 25.7% between the two surveys. Parasitaemia risk doubled in the semi-arid and low-risk seasonal zones during 2020. The greatest reductions in endemic zones were observed in Homa Bay (76.5%), Vihiga (69.1%) and Kwale (84.3%) counties. The highest increase was observed in the semi-arid seasonal zone, including Baringo (5 times), and Turkana (5–6 times).

Figure 3 presents the pixel-level relative reduction in parasitaemia risk in children aged < 5 years (A) and children aged 5–14 years (B). This figure demonstrates that there was no reduction in parasitaemia risk in the northern part of the country, the highland epidemic regions and Nairobi between the two surveys. County-specific population-adjusted parasitaemia prevalence are provided in Table A6 and Table A7 (Additional file 1) for children aged < 5 and 5–14 years, respectively.

Fig. 3

Probability of the decline in parasitaemia risk between 2015 and 2020 among children aged < 5 years (A) and 5–14 years (B)

Estimates of the number of infected children (in thousands)

This study estimated that, during the 2015 and 2020 surveys, there were 421.3 (342.2-517.9) and 296.2 (197.2-488.1) thousand infected children aged < 5 years, respectively, which translates to 29.7% reduction in the number of infected children (Table 4). The number of infected children aged 5–14 years was estimated to be 1007.0 (832.1-1182.6) and 800.8 (629.8-1075.9) thousand during the 2015 and 2020 surveys, respectively. This corresponded to a 20.5% overall reduction in the number of infected older children. The highest reduction in the number of infected children aged < 5 was observed in the low-risk zone at 50.1%, whereas for children aged 5–14 years, the highest reduction was in the coast endemic zone at 48.2%. County-specific relative reduction in number of infected children are presented in Table A6 and Table A7 (Additional file 1) for children aged < 5 and 5–14 years, respectively.

Table 4 Estimated number of infected children and population-adjusted parasitaemia prevalence by age group and epidemiological zone

Effects of intervention coverage on the changes of parasitaemia risk

Table 5 summarizes findings for the association between changes in coverage of interventions and changes in climatic/environmental factors between 2015 and 2020. An increase in nightlights was associated with a reduction in parasitaemia odds in 2020 compared to 2015 (aOR = 0.38, 95% BCI: 0.17–0.79), and (aOR = 0.66, 95% BCI: 0.41–1.04) among children aged < 5 and 5–14 years, respectively.

Table 5 Posterior estimates for the effect of interventions adjusted for difference in intervention coverage and climatic/ environmental factors between 2015 and 2020

Discussion

This study presents geostatistical analysis of the Kenya 2015 and 2020 MIS data using ZIB models to identify important intervention and climatic predictors and their net effect on malaria parasitaemia risk, produce malaria risk maps, define changes in malaria risk (including probabilities of parasitaemia risk reduction), and estimate the number of children infected with malaria parasites by age groups. The study reveals that ITN use, rainfall and nightlights indicative of urbanization and improved SES, exhibit good predictive ability for malaria parasite risk in Kenya. Overall, the study observes a reduction in malaria prevalence between the two surveys, although this reduction is notably heterogeneous, with certain areas experiencing increased parasitaemia risk, particularly the semi-arid seasonal malaria epidemiological zone. These findings offer valuable insights for the Kenyan NMCP in focussing or adjusting control efforts and aiding county health departments in implementing localized malaria control interventions. Furthermore, the ZIB model employed in this analysis holds potential for estimating prevalence in countries with low parasitaemia risk and those progressing towards elimination.

Rainfall was an important driver of parasitaemia risk in 2015, aligning with previous findings [9, 41]. An increase in rainfall was associated with increased parasitaemia risk. However, in the 2020 survey, rainfall was not an important predictor despite its increase. This could be attributed to the counter effect of increased nightlights between the two surveys. Increased nightlights indicate more urbanization and improved SES which could suggest better housing that minimizes human-vector contact. A similar analysis found no association between rainfall and parasitaemia risk in Burkina Faso, which the authors attributed to the low transmission season when the survey was conducted [42]. Though climate change may alter seasonal patterns of rainfall in the future, rainfall seasons are still within the expected months in Kenya. If implemented, and aligned to the rainy seasons, seasonal malaria chemoprevention (SMC), especially among older children who act as reservoirs, may help in reducing the burden of malaria [24].

Land Surface Temperature (LST) was not associated with malaria risk in this study. Optimal temperature is necessary for the development of malaria vectors and transmission [25]. Given the diverse nature and range – from the vast hot semi-arid north, to the highlands – the effect of temperature may have been masked. Previous studies have associated temperature with malaria risk and shown that extremely cold and hot areas may not be suitable for malaria transmission [14, 17, 41]. Changes in temperature may be associated with increase in malaria risk observed in 2020 in some of the highland counties in Kenya, including Kericho, Nyamira and Bomet. Though generally malaria parasitaemia risk is lower in higher altitudes as observed in the current and other studies [13, 17, 43].

In Kenya, like in other countries in the SSA, ITNs remain the primary malaria prevention tool [27]. In this study, we found increase in the proportion of individuals using ITNS to offer a protective effect against malaria, especially in 2020. However, an increase in ACT use was associated with an increased risk of parasitaemia. This positive effect could be due to the higher concentration of ACTs in the endemic areas compared to other areas. Our findings are not unique to Kenya; Giardina et al. observed varied effects in Angola, Liberia, Mozambique, Rwanda, Senegal and Tanzania [13]. However, ITN and ACT use were found to significantly reduce malaria risk in Uganda [41]. These findings indicate that control interventions are useful in reducing parasitaemia risk in areas where they are deployed, and with the possible expansion of malaria risk, the NMCP should consider rolling out control efforts in the areas with increased parasitaemia risk.

Overall, this study observed a reduction in parasitaemia risk between 2015 and 2020. The reduction was greater among children aged < 5 years (31%) compared to those aged 5–14 years old (26%). This study attributes the overall decline to an increase of nightlights, which is a strong indicator of improved SES and urbanization, as well as increased ITN use. It is important to note that although, in general increased nightlights are associated with reduced malaria risk, this increase may not be homogenous. This was observed in the case of Nairobi, an urban setting in a low malaria transmission zone, where there was an increase in malaria risk between the two surveys. This highlights the need for context-specific analysis to better understand localized malaria transmission dynamics. Similar findings have been observed in other studies [21, 41], where urbanization and improved SES were associated with reduction in malaria risk, which could be beneficial in Kenya’s journey towards malaria elimination. However, caution is necessary with the detection of Anopheles stephensi in cities across the Horn of Africa [44, 45], as this invasive vector is known to transmit malaria in urban areas too.

Maps depicting the probabilities of reducing parasitaemia risk and county-level estimates of prevalence reduction could assist the Division of National Malaria Program (DNMP) and county-health departments in prioritising local control measures, targeting interventions and optimising resource allocation. Of note is the significant increase in malaria risk in semi-arid seasonal epidemiological zones, especially in Turkana County, and the decline in Homa Bay County. The increase in malaria risk in Turkana County, located in north-western Kenya, may be attributed to changes in land use, which includes resource extraction and irrigation schemes that create breeding grounds for disease vectors [46] and climate suitability. Conversely, the decline in malaria risk in Homa Bay County is attributable to the IRS intervention that has been implemented in the county since 2017 [27]. The increase in malaria risk in Nairobi may be due to the importation of cases, especially from Western Kenya [47] and rapid urbanization that has led to settlement in riverine areas. Additionally, improved infrastructure over the past few years has increased movement between the regions. These findings highlight the need for intensified surveillance, deployment of diagnosis tools in expanded regions, and consideration of IRS use in high-burden areas.

With the increased pressure on health systems in SSA due to the double burden posed by infectious and lifestyle diseases, the gains in reduction of malaria may be reversed and the timelines for various milestones delayed. Additional malaria prevention and control methods may be helpful in accelerating the countries strategies towards elimination. Such interventions may include the incorporation of novel tools such as spatial repellents and next generation insecticide products such as attractive targeted sugar baits (ATSB) in the vector control toolset, depending on where more control is required [48, 49]. Adoption of the “zero malaria starts with me” campaign, an educative campaign proposed by the African Union which seeks to have malaria prevention personalized by the populations could also be beneficial in the fight against malaria [7].

This study had a few limitations. Firstly, the MIS surveys are powered for estimation of national and regional prevalence, so the raw data may not fully capture localized contexts. However, the Bayesian geostatistical models use covariates to predict prevalence in unsampled locations, thereby estimating local variation. Secondly, the surveys were conducted during different seasons, with the 2015 survey conducted during the long rainy season and the 2020 survey during the short rainy season. This seasonal misalignment may bias our findings. However, since risk factors were averaged over a period of six months prior to the end of the survey, we believe that the estimates provided here are robust. Implementing year-long rolling MIS would help address this challenge, but it remains infeasible due to high costs. Lastly, there were differences in interventions between the two surveys. In the case of the 2015 MIS clusters, they were combined with the 2014 DHS clusters to improve the precision of estimates of ITN coverage indicators. Unfortunately, this was not possible for the 2020 MIS, though the survey included more locations compared to the 2015 MIS. Nevertheless, these estimates were aggregated at the county-level, making them crudely comparable.

Conclussion

In this analysis we employed robust Bayesian geostatistical methods to assess malaria risk changes in Kenya between 2015 and 2020. The rigorous model and variable selection methods, and the modelling approach can be used in other settings with low prevalence, or when modelling sparse survey data. Given that rainfall and nightlights are important drivers of malaria risk, it is essential to align control tools seasonally and maintain vigilance, especially with urbanization. While the study observed reduction in malaria risk, some regions, particularly in the semi-arid and low-risk areas in the northern part of the country and sections of the highland areas, experienced significant increase in risk. The high spatial resolution maps depicting parasitaemia risk and probabilities of risk reduction provide valuable information for the DNMP and county government health departments. To capture spatio-temporal variations, we recommend sub-national analyses of monthly incidence data. These analyses would enable real-time monitoring of spatio-temporal changes in the distribution of malaria. This is particularly useful in advising the expansion or revision of parasite and vector control tools, and in targeting resource deployment.