Modelling of the Abundance of Malaria Mosquitoes Using Poisson Mixed Model ()
1. Introduction
Mosquitoes are insects that can cause morbidity and mortality through the diseases they can carry. Beck-Johnson et al. [1] point out that mosquitoes not only carry diseases that affect humans but are also very efficient vectors of these diseases. Mosquitoes are responsible for transmitting some of the most devastating diseases today. For instance, they may cause several problems affecting activities that have a massive economic impact such as loss of economic productivity in the workforce through mortality of staff and also pose a threat to livestock [2]. Consequently, the cumulative effect in the long term may cause a decrease in national economic capacity and development. According to the World Health Organization (WHO) estimates, the number of cases and related deaths recorded in 2017 were 219 million and 435,000, respectively [3].
Mosquitoes feed on a range of different host vertebrates and some species have developed a characteristic host preference, feeding preferentially from humans, mammals other than human and birds [4]. In Africa, A. gambiae and A. funestus are found together in many areas of the continent [5]. Mboera et al. [2] observe that in Sub-Saharan Africa, there are variations in anopheline mosquito composition and malaria transmission which are within districts and between seasons. But, Smith and others cited in [3] argue that the differences in micro-ecological and socio-economic factors included vector density heterogeneity, mosquito survival, vector host contact and their innate feeding preference are likely to have contributed to these variations. The survival of Anopheles mosquitoes is dependent on climatic conditions such as rainfall pattern, temperature and humidity.
Understanding the biology and behavior of Anopheles mosquitoes can help understand how malaria is transmitted [3]. The abundance of Anopheles mosquitoes poses a threat to the population since they are vectors that transmit some deadly diseases of public health concern, of which one notorious example is malaria. Thus, knowing where they rest is important to assess the impact of existing control methods and develop new methods to eliminate them. Hence, understanding the distribution of Anopheles mosquitoes, how they interact or differ could give a guide for a critical aspect of our future ability to control malaria. Research in this area will provide the critical scientific base for the development of strategic planning in malaria control programmes and ensure that the most effective vector control intervention tools are used.
Kreppel et al. [6] used generalized linear mixed effects models (GLMM) in investigating the variation in proportion of malaria vectors. They considered A. gambiae and A. funestus within the total catch, fitted trap type as fixed effect and night, date and house as random effects. They found that A. gambiae (which they took as Anopheles arabiensis) was caught resting mainly outdoors while A. funestus was primarily caught resting inside.
In the study of Mayagaya et al. [7], the daily abundance variations of mosquito vectors caught using different traps and resting collections, were analysed using generalized linear mixed models. Two separate analyses were carried out for A. gambiae and A. funestus, using the presence of livestock as a covariate. They found that the rate of A. funestus and A. gambiae resting indoor was lower and no significant association with the presence of livestock. They found that the proximity of livestock was associated with the abundance of A. gambiae but not with A. funestus.
Kabula et al. [4] argue that little is known about the distribution of Anopheles mosquitoes that are vectors of the disease that kills millions of people, how the species interact and differ across the country. Therefore, the abundance of the malaria vector (mosquito) would be a key determining factor for malaria risk which is a life-threatening disease. Thus, knowing their resting behavior is important for implementing control methods. The aim of this study was to investigate and determine if there was a change over time in the resting behaviour of A. gambiae and A. funestus and examine the effects of environmental variables. The remaining parts of this paper are organized as follows. We discuss the methodology considered in Section 2. The results and discussion are contained in Section 3 and lastly, we conclude the paper in Section 4.
2. Modelling Count Data
Count data refers to observations that can potentially take any non-negative integer value [8]. The aim of count data models is to explain the number of occurrences of an event. In this section, we review some of the most widely used statistical models for count data.
2.1. Poisson Model
The Poisson model can be considered as the standard model for count data [9]. The probability function of a Poisson variable
is given by,
,
y is an integer, (1)
where
and
.
Hence, the Poisson distribution is characterized by a linear relationship between its variance and mean. To model
using explanatory variables
, the log link function is used, although other choices are also available. More specifically, we write
, (2)
where;
•
is the model intercept.
•
are the regression coefficients which regulate the effects of the explanatory variables
.
•
computes the natural logarithm.
If the variance of the data is greater than the mean, the Poisson assumption is questionable. We refer to this case as over-dispersion. Ignoring over-dispersion can lead to invalid inferences on the regression relationship [10]. Probabilistic models that account for over-dispersion are the following.
2.2. Quasi Poisson Model
These models are only defined by specification of their mean, variance and a dispersion parameter k [4], that is;
and
.
We point out that a distributional form for the variable
is not required for fitting of the model and making inference on
. If k = 1, we recover the identity relationship between the mean and the variance of the Poisson model. Another advantage of the quasi Poisson model is it can be fitted using the same algorithms for the Poisson model [11].
2.3. Negative Binomial Model
The negative binomial model is obtained by assuming that the mean
of the Poisson model is Gamma distributed [8]. The density function of the negative binomial is given by,
(3)
where k is the dispersion parameter. The symbol
denotes the gamma function and is defined as
. The mean and variance of
are now
and
. The Poisson model is recovered as a limiting case for
.
2.4. Poisson Mixed Model
This model is an extension of the standard Poisson model by the inclusion of the random effects in the linear predictor. One example is given by the random intercept model which is obtained as follows. Let
for
denote random variables, having joint multivariate Gaussian distribution with mean 0 and covariance matrix
. We then assume that
, conditionally on
are mutually independent Poisson variables with mean
such that,
, (4)
It can be shown that this model has larger variance than the standard Poisson model. This model can also be further extended by allowing the regression coefficients
to vary across the observations. In Equation (4),
would then be replaced by
(assumed to follow a Guassian distribution with some appropriate mean and variance). More importantly, this model assumes that the effect of
on the response is not homogeneous but might depend on both measured and unmeasured variables.
The main difference of the outlined modelling approach with respect to the previous models, is that over-dispersion is accounted for by specification of the conditional mean of
for given latent variables (or random effects)
and/or
. Quasi-poisson and negative binomial models, instead, allow for over-dispersion through the parameter k which regulate the relationship between the marginal mean and variance of
. These two different approaches are then used to answer substantially different research questions. The former approach is useful when it is of interest to model the different sources of heterogeneity inherent to the data. The latter techniques are instead more appropriate when inference should be drawn on the general population. For these reasons, we consider Poisson mixed models to be more suitable to address the specific objectives of the present analysis.
3. Methodology
The collection of the malaria vector, A. funestus and A. gambiae mosquitoes were collected in four villages of Kilombero River Valley by the use of indoor and outdoor traps over four years from 2012 to 2015. Each village was visited several times over the four years and each time, mosquitoes were trapped inside and outside the houses. Different houses were chosen each time a visit was made to the villages and collected mosquitoes over four (4) nights in each house. Information on temperature, humidity and saturation deficit were also collected from each visited household. To catch mosquitoes inside required a different trap than to catch them outside. Traps were set in the evening at 6 pm and emptied in the morning at 6 am. For every inside trap there are 10 outside traps.
3.1. Description of Variables
The variables of interest used in this study were determined according to the study objectives and are shown in Table 1.
3.2. Model Formulation
The collected data were subjected to verification for consistency, uniformity and accuracy with the researcher in charge. Data were then exported and analyzed in the R software environment [12]. The mixed effects model was used to quantify the association of mosquito abundance with environmental factors.
Let
and
be independently and identically distributed. Let
denote the monthly counts of mosquitoes at village
at time
. We then assume that conditionally on
and
,
are mutually independent Poisson variables with mean
where
is the number of traps used to catch the mosquitoes. Using a log link-function, we then write
(5)
(6)
where,
•
is vector of dummy variables identifying each of the four villages;
•
is a vector of environmental variables;
•
is a binary indicator of the type of trap, taking value 1 if outdoor and 0 if indoor;
•
is a random intercept which accounts for extra-Poisson variation induced by unmeasured explanatory variables;
•
is a random effect accounting for unmeasured explanatory variables which affect the log-ratio between the rate of mosquitoes outdoor and indoor;
•
corresponds to the log ratio between the rate of indoor and outdoor mosquitoes at time t.
Table 1. Shows description of variables and the codes used.
The standard error is computed by utilizing the formulation,
To understand how the resting behaviour of mosquitoes changes over time, we use Equation (6) which corresponds to the log-ratio between the rate of indoor and outdoor mosquitoes and therefore, represent our target for inference. The variable selection procedure was to start with a full model including interactions with all environmental variables and type of trap, i.e. indoor or outdoor. Then, the variables with a p-value greater than 5% (i.e. considering a 95% confidence level) are removed from the model. The process is repeated until all the variables have p-values less than 5%.
4. Results and Discussion
We now look at the total number of mosquitoes collected by type of mosquito and trap in the study area, that is, for A. funestus and A. gambiae.
In Figure 1, it is observed that there are differences in the numbers of mosquitoes caught by the type of trap for A. funestus and A. gambiae. Over the four years of mosquito collection in the study area, a total of 4696 mosquitoes were collected resting indoors of which 57% were A. funestus and 43% were A. gambiae. Similarly, a total of 12,028 mosquitoes were collected resting outdoor of which 13% was A. funestus and 87% were A. gambiae.
Figure 1. Mosquitoes collected by type of mosquito and trap in the study area.
Figure 2(a) and Figure 2(b) shows the mosquito densities of A. gambiae by village for indoor and outdoor trap. It is observed that the median is different for all the villages for both indoor and outdoor trap. The highest density of A. gambiae was by indoor trap in Minepa village. We also note that there were three reported large densities by indoor trap each within Lupiro (LUP), Sagamaganga (SAG) and Kidugalo (KID) and in Minepa (MIN) by outdoor trap. We observe that by both indoor and outdoor trap for A. gambiae, the mosquito densities seem to have a similar pattern. Kidugalo and Sagamaganga village have lower densities of A. gambiae mosquitoes as compared to Lupiro and Minepa village for both indoor and outdoor trap. This can be attributed to village specific environmental characteristics. For instance, Kidugalo and Sagamaganga are on a higher altitude than Minepa and Lupiro villages.
Figure 3(a) and Figure 3(b) shows the densities of A. funestus mosquito indoor and outdoor by village. The highest density of A. funestus was by indoor trap in Kidugalo village. All the villages had the lowest density of A. funestus mosquito approximately being zero for indoor and outdoor trap. We observe large densities of Anopheles mosquitoes in villages Kidugalo, Sagamaganga and Minepa respectively by indoor and outdoor trap. In this case, we observe different patterns by village in A. funestus mosquitoes indoor and outdoor.
For both A. funestus and A. gambiae, we observe a similar pattern in the resting behaviour outdoor while the pattern seems different indoors. In almost all the villages, A. funestus was mostly caught indoor while A. gambiae was caught primarily resting outdoor. This implies that A. funestus is likely to be more exposed to insecticide treated nets (INTs) and other chemicals applied to the walls such as indoor residual spraying (IRS). However, for A. gambiae mostly resting outdoor reduces the risk of exposure to chemicals and bites humans and animals outside.
Figure 2. (a) A plot of A. gambiae mosquito density indoor by village. (b) A plot of A. gambiae mosquito density outdoor by village.
Figure 3. (a) A plot of A. funestus mosquito density indoor by village. (b) A plot of A. funestus mosquito density outdoor by village.
In order to quantify the association between mosquito abundance and environmental factors and trap location, we carried out the analysis separately for A. funestus and A. gambiae.
4.1. Model for Anopheles funestus
The first model considered was for A. funestus. After applying the variables selection procedure described in Section 3.2, only temperature range showed a significant association. The final model is shown in Table 2.
Table 2. Output for A. funestus with environmental factors.
The coefficient of Lupiro village indicates that the mosquito density is 0.7609 times the density of A. funestus mosquitoes in Kidugalo village. The density of A. funestus mosquitoes in Minepa village is 0.6178 times the density of mosquitoes in Kidugalo village. In Sagamaganga, the A. funestus mosquito density is 0.5591 times that in Kidugalo village. 40% of the density of A. funestus mosquitoes were caught outdoor implying that 60% were caught indoor. The percentage change in the A. funestus density is by 10% for every 1˚C increase in the temperature range. The interaction effect between temperature and type of trap was significant. This suggests that the effect on the density of A. funestus resting outdoors is a change by 0.8211 compared to resting indoors for every 1˚C change in temperature range.
Temperature is an important environmental factor affecting the resting behavior of A. funestus. We observe a significant relationship between temperature and the density of A. funestus. This implies that, a higher temperature range explains the change in the resting behaviour of A. funestus. It suffices to note that temperature is dependent on several factors such as location of the village, surroundings and structure of the buildings in the area.
We now carry out diagnostic checks of the residuals.
We can observe in Figure 4 that there is a fairly random pattern in a plot of residuals against the fitted values for the model of A. funestus. The random pattern indicates that the residuals appear to show homogeneity, normality and independence. This is an indication that the model provides a adequate fit to the data. As can be seen from the correlogram plot, we observe that there is no strong serial correlation in the residuals and the autocorrelation at all the lags is nearly close to zero.
4.2. Model for Anopheles gambiae
We now look at the model for A. gambiae. In the final model, temperature range, humidity range and saturation deficit showed a significant association with mosquito density. Also, there were interactions between humidity range and saturation deficit with type of trap. The estimates are given in Table 3.
Figure 4. A plot and correlogram of residuals of the model for A. funestus.
Table 3. Output for A. gambiae with environmental factors.
The coefficient of Lupiro village indicates that the mosquito density is 3.6205 times the density of A. gambiae mosquitoes in Kidugalo village. The density of A. gambiae mosquitoes in Minepa village is 5.5301 times the density of mosquitoes in Kidugalo village. In Sagamaganga, the A. gambiae mosquito density is 1.6519 times that in Kidugalo village. The A. gambiae mosquito density decreases by 0.4315 for every 1˚C increase in the temperature range. The change in the density of A. gambiae mosquitoes is 1.1444 for a 1 unit change in the average relative humidity range. A one unit change in the average monthly saturation density changes the density of A. gambiae mosquito by 17.01. The effects of the density of A. gambiae mosquitoes resting outdoors is 1.0698 times the density resting indoors for a unit change in the average relative humidity range. The effect on the density of A. gambiae mosquito resting outdoor is a change in the density by 0.3891 compared to resting indoors for a unit change in the average monthly saturation deficit.
All environment variables were statistically significant in explaining the suitability of the environment in the resting behaviour of A. gambiae. An increase in the average relative humidity implies more A. gambiae mosquitoes resting outdoor. [13] Argue that when relative humidity rises, mosquitoes survive better and is higher at night. Also, a high temperature decreases the density of A. gambiae resting indoors. [14] In their study state that A. gambiae mosquitoes have the ability to avoid high temperatures thus, move to more humid locations that facilitate their survival. Saturation deficit plays an important role in the resting behaviour of A. gambiae as it is a significant variable. The higher the average saturation deficit in the study area, the more the A. gambiae mosquitoes rest outdoor and this can also be attributed to the location of the study area.
We now check for model fit by carrying out a diagnostic plot of residuals.
As can be seen from the correlogram plot in Figure 5, we observe that there seems to be no serial correlation in the residuals and the autocorrelation at all the lags is nearly close to zero. We can also observe that there is a fairly random pattern in a plot of residuals for the model of Anopheles gambiae. The random pattern indicates that the residuals appear to show homogeneity, normality and independence and is an indication that the model provides an adequate fit to the data.
We now compare the resting behavior of A. funestus and A. gambiae by village and use of Equation (6) for the log ratio between the rate of indoor and outdoor mosquitoes. The analysis is based on a comparison of the abundance of A. funestus and A. gambiae trapped indoor and outdoor. When the value of the log ratio is zero i.e.
, there is no difference in the log ratio between the mosquitoes resting indoors and outdoors. When
implies that the log ratio of mosquito density outdoor is greater than the density indoor and conversely when
.
We start by looking at the log ratio rate between indoor and outdoor resting of A. funestus mosquitoes. The log ratio relationship is shown in Figure 6.
As can be seen in Figure 6, the patterns are different over the villages implying that there is spatial heterogeneity in the log ratio between the rate of A. funestus mosquitoes indoor and outdoor. This implies that the densities of A. funestus are different in all the villages and we observe that the mosquitoes were not collected in 2015 for some villages. In Kidugalo village, we observe that the log ratio is negative implying that the mosquito density of A. funestus was more abundant indoor than outdoor. We observe that in both Kidugalo and Minepa villages, there seems to be no clear decrease or increase in the density of A. funestus mosquito in relation to the log ratio of the resting behaviour over time. Similarly, in Sagamaganga and Lupiro villages the density of A. funestus mosquitoes was more abundant indoor. We observe that in the two villages, the fluctuations in the mosquito densities for A. funestus are seen to almost have the
Figure 5. A plot and correlogram of residuals of the model for A. gambiae.
Figure 6. Plots of the log ratio of mosquitoes resting indoor and outdoor with time for A. funestus.
Figure 7. Plots of the log ratio mosquitoes resting outdoor and indoor with time for A. gambiae.
same pattern. Also, we observe that in all villages there are no clear signs of any shift in the resting behaviour of A. funestus mosquitoes from indoor to outdoor.
Investigating the resting behavior of A. funestus mosquito over time found that it was mostly resting indoor. There were no clear variations in the resting behaviour of A. funestus in all the four villages as observed in Figure 6. These fluctuations can be attributed to malaria control strategies in the area such as indoor residential spraying, use of insecticide treated nets and other several reasons that may affect interactions with environmental variables. This indicates that despite the high usage of antibodies for malaria in the study area, mosquitoes still enter houses as evident from the abundance of A. funestus indoor.
Now, we consider the log ratio rate of A. gambiae mosquito resting indoors and outdoors. The plot is shown in Figure 7
Similarly, it can be seen from Figure 7 that the patterns over time for all the villages differ implying spatial heterogeneity in the log ratio between the rate of A. gambiae mosquitoes indoor and outdoor. We observe both positive and negative values of the log ratio between the rate of mosquitoes indoor and outdoor for all villages. The straight line represents a case when
implying that the density of mosquitoes indoor and outdoor is equal. In Lupiro and Minepa villages, we observe a similar pattern in the resting behavior of A. gambiae mosquitoes as the density of mosquitoes is more outdoor than indoor. This indicates that maybe some A. gambiae mosquitoes were resting indoors as we have observed regardless of the presence of vector control indoors. Overall, A. gambiae mosquitoes seem to be more abundant outdoor in these villages. We observe significant temporal variation in the resting behaviour of A. gambiae in all villages.
The density of A. gambiae mosquitoes in all the villages was mostly more abundant outdoor over time. This can be attributed to the fact that A. gambiae breeds in places such as rice fields, puddles, sunlit pools and bites both human and animals. For instance, the shift in the resting behavior in some of the villages can be explained by malaria control programs. In Minepa and Lupiro, the main activities are cultivation of rice and thus, A. gambiae is mainly breeding in such fields. While in Sagamaganga and Kidugalo, they have a lot of cattle and therefore this explains the shift in the resting behavior of Anopheles to biting animals outdoor.
5. Conclusion
We conclude that A. gambiae and A. funestus show a remarkable temporal variation in their resting behaviour. A. funestus was the most abundant malaria vector indoor while A. gambiae was more abundant outdoor in the study area. The finding is in agreement with the results of [6]. Over four years of mosquito collection, 57% of mosquitoes collected indoor were A. funestus while 43% were A. gambiae. Similarly, 87% of the mosquitoes collected outdoor were A. gambiae and 13% were A. funestus. Minepa and Lupiro villages had large densities for both A. gambiae and A. funestus indoor and outdoor as compared to Kidugalo and Sagamaganga village. This can be due to their location as they are on the lower altitude of the Valley compared to Kidugalo and Sagamaganga, with possible presence of waterlogged areas. Also, the prominent activity in Minepa and Lupiro villages is farming while keeping cattle is the prominent activity in Kidugalo and Sagamaganga villages. Thus, the rice fields may be breeding sites for mosquitoes. There were variations in the resting behaviour of A. gambiae and A. funestus over time in all the villages. The patterns in the log ratio between indoor and outdoor densities of A. gambiae and A. funestus were different implying spatial heterogeneity in their resting behaviour. A. funestus was found to rest mainly indoor over time in all the villages. We found that A. gambiae was more abundant outdoor although there were some resting indoors in all the villages. Environmental variables had an impact on the abundance of A. gambiae and A. funestus. Temperature was found to significantly affect the resting behaviour of A. funestus. Whilst humidity range, temperature and saturation deficit had an effect on the resting behavior of A. gambiae mosquitoes in the study area. The findings of this study could be used as a guide to implementing appropriate intervention measures. The use of indoor residential praying (IRS), insecticide treated nets (ITNs) and other insecticides may then impact on the abundance of A. funestus. However, for mosquitoes outdoor, the use of mosquito repellents and other interventions would help reduce the abundance of A. gambiae which was more abundant outdoor. Ensuring that all possible mosquito breeding sites are eliminated is also a way to eliminate mosquitoes outdoor as this is prevention through source reduction.
Acknowledgements
The authors would like to express their gratitude to Ifakara Health Institute for providing the data.