The most common indicator to evaluate vector control interventions such as LLIN and IRS relies on malaria transmission through estimation of the Entomological Inoculation Rate (EIR). EIR is the product of the Human Biting Rate (HBR; number of bites of malaria vectors per human per unit time) and the prevalence of Plasmodium infection in mosquitoes. HBR is usually measured using the Human Landing Catches (HLC) counting technique that is the method of reference to quantify the human-vector contact [3].

In 28 villages of Southern Benin, a recent cluster randomized controlled trial (RCT) aiming at comparing the efficacy of combined LLIN and carbamate IRS or carbamate-treated plastic sheeting (CTPS) with a background of LLIN coverage did not show benefits of the combination for reducing HBR and EIR [4]. In the study area, high variations in the density of malaria vectors were observed in time and space [5] and there were many localities with zero mosquitoes collected during several nights.

The most ancient and popular statistical distribution used to describe count data is the Poisson distribution that assumes equidispersion of the counts. However, in real datasets, these counts are often overdispersed [6]–[8] and there are various means to demonstrate it [9]–[11]. Among the causes of overdispersion is the excess of zeros. Within the context of malaria vectors counts, the excess of zeros may result from the absence of mosquitoes at some locations (houses, village…) or during some period of time (dry season, cold temperatures…).

To deal with such overdispersed data with excess zeros, Johnson and Kotz [12] introduced the zero-inflated Poisson model (ZIP); i.e., a Poisson mixture model that combines a point mass at zero with a Poisson count distribution. Later, Lambert [13] extended this model to allow for covariates. Another way to deal with count overdispersion is the use of the negative binomial (NB) model or, better, the zero-inflated negative binomial (ZINB) model constructed on the same principle as that of the ZIP.

Besides these well-known models, other finite mixture distribution models have been proposed (e.g., McLachlan and Peel [14]) and have been the object of numerous applications. In fact, these models extend the previous ones; instead of considering a mixture of two distributions as with the ZIP or the ZINB, they consider a mixture of three or more Poisson or NB distributions. In addition, a non-parametric approach of the maximum likelihood introduced by Aitkin [15] has shown to be an excellent tool to allow for overdispersion. An extension of this approach by the same author [16] allowed its application to repeated measurements. Thus, a non-parametric mixture of Poisson model (NPMP) seems adapted to take into account the frequent changes in vector counts in various sites of a study zone.

McCullagh and Nelder [29] asserted that whenever the variable of interest is a count, its distribution is often an overdispersed Poisson distribution. The present data are another illustration of this assertion. The first part of this work aimed at comparing which distribution among Poisson, NB, ZIP, ZINB and NPMP better fit on counts of malaria vectors recorded using the HLC technique. Both NB and NPMP models dealt with the excess of zero, with overdispersion and provided the best predictions of the distribution of the observed data. However, unlike NB model, the NPMP does not do any further assumption about the distribution of the means of malaria vectors counts. Besides, the hierarchical structure of the observed data was taken into account by a NPMP model conditional on “village”. Based on a posterior probability criterion, the NPMP model allowed ranking the villages in four latent classes according to the mean of vector density after adjustment on environmental and climatic covariates. The optimal number of latent classes was established on conventional criteria. Furthermore, the part each covariate played in the variability of malaria vector density in the area was estimated by the fixed effects of the model. However, the present study could not take into account all the possible hierarchical levels of the data because of the limits of software R in dealing with latent classes. Indeed, function allvc of package “npmlreg” cannot deal with more than two levels. We considered thus the catches at all sites of the same village as repeated measurements of the same variable. Therefore, we were not able to take into account the possible correlation between the counts from houses within the same village [30]. “Human bait” is another level that could induce correlation in the data but there is no sufficient information about all mosquito collectors. Besides, the rotation of the collectors during data collection reduces considerably such a correlation. “Season” could be another possible level of correlation; it was taken into account through rainfall data which is the main seasonal factor in the context of malaria vector density.

Moreover, the numbers of collected vectors during the 8 surveys are assumed to be uncorrelated although one may speculate about a correlation structure along time. Nevertheless, the correlation between mosquito counts from successive surveys is deemed to be very low because the time span between two successive surveys is 6 weeks whereas the lifespan of the vectors is only 3 to 4 weeks. Studying the correlation between counts from two nights during the same survey may reveal interesting results.

In southern Benin, both spatial and temporal heterogeneities in vector densities were mentioned by Djènontin et al. [5]. This can be explained by some factors we found associated with the density of malaria vectors. Firstly, cumulated rainfalls during the 15 days preceding the catches were positively associated with vector density as previously reported in Benin [30]. Moreover, the mean number (over all surveys) of rainy days was positively associated with the vector density whereas the deviation at each survey from this mean was negatively associated with the vector density. This suggested that high frequency of rainy events might flush out vectors breeding sites [31]. The vector density was lower in villages with water supply; this could be due to the absence of water storages that could have provided breeding sites for malaria vectors [32], [33]. Moreover, the presence of irrigated market gardening could have provided breeding sites [34], [35] and then, increased the density of vectors in villages closed to this activity as previously observed in Benin [36]. Permanent freshwaters of the Toho Lake could also have provided breeding sites for both An. funestus and An. gambiae [37]–[39] that are both present in our study area [5]. This explains why the vector density decreased when moving away from freshwater bodies as showed by Amek et al. [40] in Western Kenya. The presence of cattle was negatively correlated with vector density suggested that a part of the vector population could have bite on cattle instead of human. More vectors were caught in multi-cluster villages than in single-cluster villages. This might indicate that a multi-cluster village layout might increase the attractiveness of the village for malaria vectors because of the extra vegetation surrounding houses. Thus, the attractiveness of a multi-cluster village may be higher than that of a single-cluster village of same size. Catches were also more abundant outside than inside the houses. This indicates an exophagic behavior of malaria vectors in the study area. As suggested by two studies in the OKT region [4], [41], a part of the exophagic population of vectors could have avoided indoor residual insecticides.

One unexpected finding of the present study was that the NDVI was negatively correlated with the density of malaria vectors. This finding contrasts with several studies that used satellite imagery at a lower resolution [42], [43] but agrees with a study carried out in Burkina Faso that used the same SPOT images than ours [44]. In this study, the authors found a negative relationship between the larval productivity in ponds and the NDVI calculated from high resolution SPOT images. Indeed, a high NDVI might reflect the presence of submerged vegetation or water covered with vegetation that are usually related to very high Anopheles larval densities [44]–[46]. Moreover, the NDVI usually decreases with freshwater and unvegetated surfaces likely to provide breeding sites for the malaria vectors [37], [47]. Nevertheless, the discussion about the NDVI effect can be more complex because of the co-existence in the region of two major malaria vectors with different breeding-site requirements.

In this work, villages were ranked into four classes of increasing mean malaria vector density but we were not able to find any relationship between this grouping structure and the vector control intervention implemented in the village. This confirms the finding of Corbel et al. [4] who demonstrated with the same data, that vector density was not significantly different between treatment arms (TLLIN, ULLIN, TLLIN+IRS, and ULLIN+CTPS).