The relative contribution of climate variability and vector control coverage to changes in malaria parasite prevalence in Zambia 2006–2012

Four malaria indicator surveys (MIS) were conducted in Zambia between 2006 and 2012 to evaluate malaria control scale-up. Nationally, coverage of insecticide-treated nets (ITNs) and indoor residual spraying (IRS) increased over this period, while parasite prevalence in children 1–59 months decreased dramatically between 2006 and 2008, but then increased from 2008 to 2010. We assessed the relative effects of vector control coverage and climate variability on malaria parasite prevalence over this period. Nationally-representative MISs were conducted in April-June of 2006, 2008, 2010 and 2012 to collect household-level information on malaria control interventions such as IRS, ITN ownership and use, and child parasite prevalence by microscopic examination of blood smears. We fitted Bayesian geostatistical models to assess the association between IRS and ITN coverage and climate variability and malaria parasite prevalence. We created predictions of the spatial distribution of malaria prevalence at each time point and compared results of varying IRS, ITN, and climate inputs to assess their relative contributions to changes in prevalence. Nationally, the proportion of households owning an ITN increased from 37.8 % in 2006 to 64.3 % in 2010 and 68.1 % in 2012, with substantial heterogeneity sub-nationally. The population-adjusted predicted child malaria parasite prevalence decreased from 19.6 % in 2006 to 10.4 % in 2008, but rose to 15.3 % in 2010 and 13.5 % in 2012. We estimated that the majority of this prevalence increase at the national level between 2008 and 2010 was due to climate effects on transmission, although there was substantial heterogeneity at the provincial level in the relative contribution of changing climate and ITN availability. We predict that if climate factors preceding the 2010 survey were the same as in 2008, the population-adjusted prevalence would have fallen to 9.9 % nationally. These results suggest that a combination of climate factors and reduced intervention coverage in parts of the country contributed to both the reduction and rebound in malaria parasite prevalence. Unusual rainfall patterns, perhaps related to moderate El Niño conditions, may have contributed to this variation. Zambia has demonstrated considerable success in scaling up vector control. This analysis highlights the importance of accounting for climate variability when using cross-sectional data for evaluation of malaria control efforts.


Background
Scale-up of vector control interventions, specifically ownership and use of insecticide-treated nets (ITNs) and indoor residual spraying (IRS), has been shown in multiple settings to reduce malaria morbidity and mortality [1][2][3][4][5]. Zambia made significant progress in scaling up ITN and IRS coverage from 2005 through 2012, with progress documented through repeated nationally-representative household surveys [6,7]; over 10 million ITNs were distributed between 2005 and 2012, and roughly 850,000 structures were sprayed annually during this period. Nationally-representative Malaria Indicator Surveys (MIS) were conducted in 2006, 2008, 2010 and 2012 to monitor the progress in household coverage of ITNs and IRS at the population level, as well as to provide estimates of malaria parasite prevalence in children 1-59 months at these time-points.
Although household ownership and use of ITNs and IRS coverage as measured by these surveys increased nationally from 2006 through 2012, the proportion of children < 5 years old with a malaria parasite infection declined from 22.4 % in 2006 to 9.3 % in 2008 [8], but increased somewhat to 15.9 % in 2010 and 14.6 % in 2012. The increase between 2008 and 2010 was variable across the country and greatest in rural areas, nearly doubling in some provinces; the proportion of children with severe anemia increased in a similar pattern with parasitemia between 2008 and 2010, from 4.3 to 9.2 %.
Identifying possible causes of this observed resurgence in infection prevalence between 2008 and 2010 has broad implications for malaria control programs across Africa. A historical analysis of global occurrences of malaria resurgence suggests that in most cases a flagging of malaria control effort led to an increase in incidence rates [9]. There is growing evidence of potential ITN effectiveness decay in sub-Saharan Africa (SSA) due to insecticide resistance and changes in vector biting behavior [10][11][12]. There is further speculation that anomalous weather patterns-perhaps related to a relatively strong El Niño episode in late 2009 and early 2010-may have influenced this resurgence in transmission, as observed elsewhere [13,14].
Weather variability is a well-known driver of malaria transmission [15]. High-resolution climate and environmental data are being used with increasing sophistication in geostatistical modeling frameworks for malaria risk mapping purposes [16], where the objective is usually the production of mean endemicity surfaces [17], incidence prediction [18] or to examine changes in malaria parasite prevalence over time [19,20]. However, until recently climate data have seldom been directly incorporated into evaluations of program impact [21][22][23]. Three recent examples where climate data were successfully incorporated include the evaluation of vector control scale-up and incidence data in Eritrea by Graves and colleagues [21], work by Giardina and colleagues comparing changing parasite prevalence and vector control coverage in five countries [24], and a continental-scale evaluation by Bhatt and colleagues [5]. In order to obtain unbiased estimates of the impact of malaria control in program evaluations and assess any potential change in their effectiveness over time, analyses must directly incorporate changes in climate that influence malaria transmission potential.
Zambia reported increases in malaria parasite prevalence and health system reported clinical incidence between 2008 and 2010 that continued in 2012, despite continued scaleup of malaria control interventions during this time. This paper assessed the association between inter-annual climatic and other environmental variables, IRS and ITN coverage, and changes in malaria parasite prevalence between 2006, 2008, 2010 and 2012 while accounting for confounding factors at the subnational level. We used geostatistical models to estimate the relative contribution of subnational changes in IRS and ITN coverage and climate to changes in malaria parasite prevalence over this period.

Malaria indicator surveys
We used data from three Malaria Indicator Surveys (MISs), each of which was conducted at the end of the high transmission season between April and June in 2006, 2008, 2010 and 2012. The sampling design and questionnaire for these surveys has been described elsewhere [25]. Briefly, the sample size and standard enumeration areas (SEAs) were selected to provide precise estimates of ITN coverage at the national, provincial, and urban/rural levels. At the time of the surveys, there were 9 provinces and 72 districts in Zambia. While a new province was demarcated in 2011, we maintained the original 9 provinces for consistency in this analysis. A two-stage sampling design was used, with the primary sampling units consisting of standard enumeration areas (SEAs) selected proportional to the estimated population size (PPS) of each within provincial and urban/rural strata. Within each selected SEA, field workers conducted a complete household enumeration using personalizeddigital assistants (PDAs) equipped with GPS, and selected 25 households for questionnaire administration to the household head and caregivers of children under 5 years of age. Latitude and longitude were collected for each household. For SEA-level geographic information, we determined the centroid for each SEA by averaging household latitudes and longitudes.

Primary outcome
Malaria parasite prevalence in children 1-59 months served as the primary outcome. Malaria parasite prevalence was ascertained from infection status by quality controlled slide microscopy for all children 1-59 months in selected households during each MIS; HRP2 rapid diagnostic tests (RDTs) (ICT Malaria Pf ) were used in the field to provide point of care diagnosis and treatment for children testing positive.

Measures of explanatory variables and potential confounding factors
Publicly-available remote-sensing climate data referenced to the centroid of each selected survey SEA were compiled from the following sources: the African Data Dissemination Service for 10-day rainfall estimates at 8 km resolution [26], the USGS Hydro 1 k dataset for elevation at 1 km resolution [27], the MODIS satellite data repository for 16-day enhanced vegetation index (EVI) at 1 km resolution [28], the monthly Temperature Suitability Index (TSI) from the Malaria Atlas Project (MAP) [29], and Worldpop for population estimates at 1 km resolution [30]. We also investigated program administrative shape-files (.shp) for road networks, health facilities, and water bodies. Time-variant climate and environmental data (rainfall, TSI, and EVI) were extracted for the four months preceding the start of each survey at 1 km resolution.
In each MIS, a net roster was used to count the number of mosquito nets present in each selected household, including when the net was procured and who in the household slept under the net the previous night. If and when IRS was conducted was ascertained based on the recall of the head of household during each MIS. We defined a mosquito net as an ITN if it was a long-lasting insecticide-treated net (LLIN), or any net that had been insecticide treated within the past 12 months [31]. Both household and community-level ITN ownership were considered: for community-level coverage, we calculated the ITN density (number of ITNs per household members) for each SEA. We defined IRS as the survey respondent reporting that household walls were sprayed within the 12 months preceding the survey date. Similarly, we calculated the proportion of households having received IRS per SEA as a measure of community-level coverage.
As a measure of household socioeconomic status, principle-components analysis (PCA) was used to create a wealth index from a list of household assets captured during each survey [32]; this asset list did not include any household vector control items. Population density was calculated from the Worldpop raster and used as a measure of urbanness. We classified areas with < 1,000 population per km 2 as rural areas and areas with > 1,000 population per km 2 as urban/peri-urban areas. This has been shown to be a reliable cut-off for stratifying malaria risk [33]. We also calculated the Euclidean distance from each household to the nearest water body, including permanent lakes and major rivers.

Variable selection and statistical analysis
Descriptive comparisons of household vector control and parasite prevalence were conducted across surveys both nationally and by province, using survey-weighted point estimates.

Climate variables
For modeling and predicting parasite prevalence using time-varying climatic data, we considered both biological plausibility and statistical fit to select climatic variables based on lag periods preceding household data collection. For EVI, lags up to the beginning of survey data collection were considered as this index represents the effect of cumulative rainfall and humidity [further details of the selection of climate variables are presented in Additional file 1]. The best-fitting lags within biologically plausible time frames included in final models were the 20-day period of rainfall ending 6 weeks prior to the start of data collection, and the 16-day EVI in the week preceding the start of survey data collection.

Model selection
Bivariate and multivariable associations between the outcome of parasite infection and explanatory variables were assessed using random effects logistic regression models for each survey round. Linear, categorical, and non-linear relationships between the set of explanatory variables and parasite prevalence were assessed in separate models for each survey; interaction terms between survey year and ITN coverage, climate, and other covariates were assessed in final pooled models. Final non-spatial model selection was based on the Akaike's Information Criterion (AIC); final spatial model selection was based on the Deviance Information Criterion (DIC). Model fit assessments were conducted in Stata 13.1 [34] and R 3.3.0 [35].
To account for spatial autocorrelation and to produce continuous prediction surfaces for examining spatial patterns in ITN coverage and parasite prevalence, final models were run in a Bayesian geostatistical framework with spatially-correlated random effects (further details in Additional file 1). We compared results of several different geostatistical models including individual, household, and SEA-level parameters to confirm that effect estimates used for final prediction models were not modified by individual level parameters. Additionally, we compared results of geostatistical models with nonspatial random effects models. In a model ('full model') including individual, household, and cluster level parameters we included individual and household-level sociodemographic and vector control parameters (ITNs and IRS), and SEA-level measures of vector control (ITNs and IRS), geography (elevation and Euclidean distance to nearest water body), and climate (rainfall and EVI). In the spatial prediction model ('prediction model') we included only the SEA-level ITN, IRS, geographic, and climate parameters. Additionally, to compare results of varying ITN and IRS coverage and climate inputs on predictions, we created separate geostatistical models of ITN coverage, defined as the ratio of ITNs:HH members, and IRS coverage for inclusion in final prediction models. All model parameters were estimated using Integrated Nested Laplace Approximation (INLA) in R [36] (Details in Additional file 1).

Model validation and prediction
We conducted validation tests on the spatial prediction model and a comparative non-spatial model by calibrating each model on an 85 % training dataset for predicting to the remaining 15 % of locations. The training and prediction datasets were selected through stratified sampling to ensure equal allocation by survey year. We calculated the mean prediction error (MPE) as a measure of bias and the root mean square error (RMSE) as a measure of accuracy. For the final prediction model surfaces, we plotted the posterior mean prediction values for each 5 km pixel (roughly 40,000 pixels covering Zambia) using INLA in R [36] and ArcGIS [37]. We evaluated several different prediction scenarios by including different combinations of climate and ITN and IRS coverage surfaces for each year as prediction model inputs. For alternative prediction scenarios we focused primarily on changes between 2008 and 2010 as this period was characterized by anomalous climate conditions and the largest changes in measured prevalence. For each prediction scenario we calculated final population-averaged parasite prevalence values per province by multiplying the mean posterior prediction raster by the Worldpop population raster for children < 5 to obtain the estimated number of infected children < 5 by province, and then dividing by the total province < 5 population. To approximate uncertainty in province level estimates, we calculated the population-averaged parasite prevalence for the 2.5 and 97.5 % posterior quantiles.

Results
The 2006 MIS covered 2,889 households, with slide results available for 1,787 children 1-59 months of age in 109 SEAs (Fig. 1). The 2008 MIS covered 4,142  (Table 1); combined, the percent of households with any vector control (ITN and/or IRS) fell significantly in both provinces. In other provinces -Eastern, Copperbelt and Central, household ITN ownership increased slightly. Between 2010 and 2012, the percent of household with an ITN and/or IRS increased substantially in Luapula (from 60.9 to 91.4 %, P < 0.001), Eastern (from 77.8 to 90.4 %, P < 0.05), and Northern (from 66.1 to 78.6 %, P < 0.05), but decreased in Western (from 79.3 to 56.0 %, P < 0.05) and Central (from 77.7 to 65.2 %, P < 0.05).
The ITN-to-population surfaces produced as inputs for spatial prediction showed that this ratio increased overall across provinces from 2006 through 2010, but with substantial spatial heterogeneity, and increased dramatically in 2012 (Additional file 1: Figure S4). The number of ITNs per population increased overall between 2006 to 2008, but most notably in the north and east of the country. Between 2008 and 2010, the number of ITNs per population remained at similar levels in most of the west and south, but fell in the north, most notably in Luapula province. In 2012, this ITN-topopulation ratio was highest in the northeast. IRS coverage was highly focal and higher near urban areas in 2006 and 2008, then increased nationally in 2010 and 2012 (Additional file 1: Figure S5 Table 2). Prevalence decreased slightly, but not significantly, in Northwestern (from 14.1 to 6.0 %, P = 0.150) and Southern Province (from 7.7 to 5.7 %). National prevalence in 2012        (Fig. 3). In 2010, resurgence of parasite prevalence was most notable in Luapula province and portions of Eastern province. From these predicted surfaces, population-adjusted parasite prevalence estimates were similar to the MIS estimates, with national and provincial level changes by survey year consistent with changes estimated by MIS tabulations (Table 4). Table 4

Discussion
We used a geostatistical model to assess the subnational associations between population-based measures of vector control coverage, inter-annual climate variability, and parasite prevalence in Zambia in 2006, 2008, 2010, and 2012 from the MIS in these years. In doing so, we were able to isolate the contributing effects of changes in ITN and IRS coverage and climate on changes in parasite prevalence at the subnational level in Zambia over this period, with a specific focus on changes that occurred between 2008 and 2010.   due to climate variability. In geostatistical prediction models, we estimated that while the national increase in populationadjusted predicted parasite prevalence between 2008 and 2010 was from 10.4 to 15.3 % (a relative increase of 47 %), prevalence would have risen to 15.7 % (a relative increase of 51 %) if vector control was the same as in 2008. We estimate that while population-adjusted parasite prevalence increased relatively by almost 80 % from 2008 to 2010 in Luapula, prevalence would have increased by close to 40 % relatively due to lower ITN coverage alone. However, given the wide uncertainty in predicted prevalences, the credibility intervals for predictions overlapped, which limits our ability to make statistical conclusions about these results. Although heterogeneity in ITN coverage over this period played a substantial role as ITN coverage fell in several of the provinces that experienced the greatest increase, these patterns do not completely explain the observed rebound. Higher community ITN density (ITNs per individual) and higher IRS household coverage were each significantly associated with reduced odds of parasite infection across all surveys, and coverage of each of these fell significantly in Northern and Luapula province between 2008 and 2010. Provinces where community vector control coverage in 2010 was at least as high in 2008 did not experience significant increases in infection prevalence, with the exception of Eastern province. In 2012, ITN coverage reached very high levels in Luapula and Eastern, but while prevalence dropped from 2010, it remained high. In Western province, ITN coverage fell in 2012, and this was associated with a more than doubling in parasite prevalence.
While IRS rates increased, the proportion of households with a newer ITN (< 2 years old) decreased dramatically in Northern, Luapula, and Eastern province between 2008 and 2010. This was likely the result of the staggered LLIN distribution campaigns that resulted in some provinces not receiving new nets until after 2010. This highlights the need for programs to work diligently at ensuring vector control coverage does not lapse in certain areas due to issues with net redistribution or spray efforts.
Climatic factors influencing malaria transmission, including rainfall, temperature suitability, and vegetation cover (EVI), were positively associated with malaria parasite prevalence in children across space and time. These changes in rainfall and temperature patterns may reflect regional effects of the El Niño Southern Oscillation (ENSO), as early 2010 was characterized by a moderately strong El Niño event, and early 2008 by a moderate La Niña event. However, the effects of the ENSO on the southern Africa region are varied, and likely dependent upon latitude [38]. Previous research found that La Niña years predict higher incidence in Botswana, which is similar ecologically to parts of southern Zambia [18]. Conversely, El Niño years have historically affected the northern part of the region [38]. This pattern may partially explain the rebound in the provinces in the northern part of country. Evidence of the influence of climate conditions favorable for malaria transmission in 2010 in southern Africa was not limited to Zambia. Similar trends were observed in subnational surveys in Malawi between 2005 and 2010 (PMI Malawi Impact Evaluation, personal communication). An increase in parasite prevalence in other countries in central Africa (e.g. Rwanda) may have also been influenced by a combination of intervention fall-back and climate factors. Trends in health facility data in Zambia and Malawi also coincide with a resurgence in malaria that began in 2009 and increased through 2010 [39].
Climate and ITN variables do not completely explain observed subnational changes in malaria parasite prevalence between 2008 and 2010. Immunological factors and other transmission dynamics between survey years may have influenced prevalence in ways not fully captured in geostatistical models. For example, older children in 2010 would have likely been less exposed to infectious mosquitoes in 2008 following the ITN distribution and a drier late-transmission season, and therefore possibly have weaker acquired immunity in 2010 than similarly aged children would have in 2008. They might therefore have been less able to clear infections or suppress parasite densities [40][41][42]. Other factors that must be considered in the resurgence in 2010 in some areas include ACT drug shortages and insecticide resistance. However, reported ACT levels were largely consistent between 2006 and 2010, and treatment seeking and treatment of fevers with an ACT has increased since 2006 [43]. Additionally, prevalence remained high in 2012, even as ACT coverage increased. As a result, these factors were not likely to have played a substantial role in increased prevalence over this period. Insecticide resistance may have contributed, but its role on a national scale is largely unknown, as studies have been very localized [10][11][12]. Mathematical modeling has furthermore suggested that the duration of the effects of LLIN provision may vary with transmission setting, such that the effective lifetime is shortened in higher transmission settings [44]. However, in our study ITNs at the community level and IRS at the household level were shown to be protective against malaria parasite infections across survey years, with no statistical evidence of ITN effectiveness decay as assessed by interaction terms between survey year and both household ITN ownership and community ITN density (non-significant interactions).
As with any cross-sectional study, this analysis was limited in its ability to adequately examine longitudinal trends. The utilization of repeated cross-sectional surveys allows only for examinations of contemporaneous association with limited causal inference [45]. Given that control with ITNs and IRS has been scaled-up nationwide, we lacked a true control group that would allow us to make stronger counterfactual arguments. We attempted to mitigate some of these limitations through the use of matching based upon residence status in non-spatial models and found associations similar to those in geostatistical models. Similarly, we approximated province level uncertainty in predictions by using the 2.5 and 97.5 % posterior quantiles at the pixel level. Given the known computational difficulties with producing aggregate-level uncertainty from geostatistical surfaces [46], we were not able in this study to incorporate the full pixel-level covariances. Additionally, further work is needed to parameterize potentially nonlinear relationships between ITN coverage and parasite prevalence in a geostatistical framework.
The use of lagged climate variables to predict prevalence is confounded by the length of a malaria infection, which averages over 200 days without treatment [47]. However, as diagnostic and treatment capacity increased dramatically over this period, prevalence would have likely tracked incidence more closely as infections were cleared with more rapid effective treatment. Similarly, as transmission drops, as it did between 2006 and 2008, some children may lose acquired immunity, and therefore present with more symptomatic disease necessitating treatment. These effects would seem to bias toward shorter infections, which could potentially allow for more accurate estimation of climatic time lags with prevalence data. A strong relationship between climate and prevalence through time is supported by recent analyses of community-level surveillance data in Eastern province showing that the seasonality in prevalence closely tracks incidence [48].

Conclusions
There is increasing effort to conduct repeated national surveys to generate data for malaria program evaluations. Without carefully considering and explicitly modeling confounding factors, especially climate, results of these surveys may lead to spurious conclusions of program effect or lack thereof, especially when assessing changes at the national level. In Zambia, it was surmised that ITNs began to fail after 2008 and contributed to increases in national infection prevalence by 2010. However, once one explicitly models changes in inter-annual climate variability at the subnational level, in conjunction with changes in ITN population coverage at the subnational level, we show that without the ITN coverage already in place in 2010, the malaria burden as measured by predicted parasite prevalence would have been worse in 2010 than actually observed. In doing so, this analysis represents a novel approach to applying geostatistical methods for malaria program evaluation. These results are important as malaria control programs attempt to interpret changes in prevalence and attribute them to program effort, or lack thereof, in the context of inter-annual climate variability. Subnational changes in malaria parasite prevalence are highly influenced by climate over short time scales and must be accounted for in assessing the effectiveness of malaria control programs.

Additional file
Additional file 1: Supplementary information: Additional modeling results and maps of vector control coverage. (DOCX 1 mb)