Ce document de travail a pour objectif de présenter une méthode de modélisation de différents traits et préférences (bioécologiques, comportementaux et phénotypiques) des principales espèces vectrices du paludisme à échelle paysagère en milieu rural ouest-africain, en vue de la rédaction d’un article scientifique portant sur ces questions. L’objectif principal de cet article est de tenter d’expliquer ces traits avec des variables explicatives diverses (climatiques, occupation du sol, de comportement humain, etc.).
Les variables que nous souhaitons expliquer sont les suivantes :
Nous avons réfléchi ce document méthodologique par rapport à la variable a expliquer “densité agressive” des vecteurs. Parmi les 4 variables mentionnées, c’est probablement la plus délicate à modéliser au regard de i) sa distribution statistique complexe (binomiale négative avec excès de zéro) et ii) le nombre de variable explicatives considérées. Nous pensons que la méthode détaillée dans ce document peut être étendue aux autres variables, qui présentent des distributions statistiques plus simples et un nombre de variables explicatives moins importantes.
La variable que l’on cherche à modéliser est la densité agressivité des espèce majeure d’Anophèles. Il s’agit donc d’une variable de comptage (valeurs numériques discrètes). Nous utilisons des données receuillies par capture sur sujet humain, sur deux zones d’études différentes - Côte d’Ivoire (CI) et Burkina Faso (BF) - pendant environ 1,5 années.
Nous proposons d’étudier séparement les deux zones d’études pour chaque variable à modéliser car i) la densité agressive est significativement plus importante en CI qu’au BF et ii) certaines variables explicatives sont disponibles dans une zone mais pas dans l’autre. Ainsi il sera intéressant de tenter de comprendre ce qui explique les différences entre les deux zones.
La distribution statistique de chacune de ces variables à modéliser est donnée dans la section “Choix du modèle statistique pour l’analyse de nos données”. L’analyse de cette distribution permettra d’orienter le choix du modèle à utiliser pour modéliser chacune des variables.
La carte suivante présente la localisation des villages de l’étude dans les deux zones. L’objectif principal du projet REACT est d’étudier, par essai randomisé contrôlé (ERC), l’impact de la mise en place de mesure de lutte anti-vectorielle (LAV) complémentaires à l’utilisation des moustiquaires imprégnées d’insecticides sur la population de vecteurs et l’incidence du paludisme. Ainsi au cours de l’étude chaque village s’est vu attribué au hasard une méthode de LAV complémentaire. Un lot de village témoin (“contrôle”) n’a pas reçu de méthode de LAV complémentaire, conformément à la méthode ERC.
A titre d’exemple, le tableau ci-après donne les dates des missions entomologiques dans les deux zones d’étude et les captures totales d’Anophèles, et les figures qui suivent donnent la distribution du nombre de piqures par homme et par nuit (human biting rate, abbrégé ma), toutes espèces d’Anophèles confondues, sur les deux zones d’études, tous traitements confondus puis en fonction des différents traitements.
codepays | nummission | date | ma |
---|---|---|---|
BF | 1 | 2017-01-23 | 282 |
BF | 2 | 2017-03-25 | 290 |
BF | 3 | 2017-05-20 | 1406 |
BF | 4 | 2017-11-09 | 458 |
BF | 5 | 2017-12-08 | 234 |
BF | 6 | 2018-01-21 | 101 |
BF | 7 | 2018-03-24 | 225 |
CI | 1 | 2016-09-30 | 34233 |
CI | 2 | 2016-11-27 | 3365 |
CI | 3 | 2017-02-23 | 969 |
CI | 4 | 2017-04-20 | 5647 |
CI | 5 | 2017-10-27 | 9817 |
CI | 6 | 2017-12-04 | 1325 |
CI | 7 | 2018-01-25 | 168 |
CI | 8 | 2018-03-24 | 2195 |
L’analyse visuelle de ces deux derniers graphiques laisse penser que la variable à expliquer (nombre de piqures d’Anophèles) suivra une distribution asymétrique, avec de nombreux comptes nuls ou proches de zéro, et une proportion moins importante de comptes positifs >> 0 . Nous vérifions cette hypothèse en représentant la fréquence des comptes (toutes missions entomologiques confondues) :
Nous confirmons donc notre hyopthèse : il s’agit de distributions binomiales négatives. Par ailleurs, nous avons vérifié que ces distributions présentent également un excès de zéros (zero inflation) en comparant la sortie de modèles GLMM ZINB et ZIP et en vérifiant que la différence entre les deux modèle était significative.
En modélisation statistique “traditionnelle”, ces types de distributions sont souvent analysées avec des modèles paramétriques type «Zero Inflated Negative Binomial » (ZINB) ou modèles de haie (hurdle model). Ces modèles sont largement reconnus par la communauté scientifique et ont été utilisés à de maintes reprises pour modéliser des compte de piqures (aggréssivité) de vecteurs (refs). Dans nos travaux, nous proposons d’utiliser une autre classe de modèles, à savoir des modèles non-paramétriques (algorithmiques). Nous justifions l’utilisation de ces modèles dans la section ‘Méthode de modélisation statistique’.
Les variables explicatives qui serviront dans les modèles mathématiques sont données dans le tableau suivant, avec les informations notamment sur les données sources utilisées et les pré-traitements effectués. Il s’agit principalement de variables climatiques (températures, précipitations, etc.), paysagères (occupation du sol), anthropiques et démographiques (population). Le tableau est aussi disponible à l’adresse suivante https://docs.google.com/spreadsheets/d/1dIeSOa2WinXvOQGLmIjA0gFdsHnb6zMMsDME-G5pyMc/edit?usp=sharing onglet “var_explication”.
Ces variables peuvent être séparées en 3 groupes :
Les variables spatialisées ont été extraites autour de chaque point de capture selon plusieurs zones tampon, afin de mesurer l’éventuelle influence de la variabilité spatiale de ces variables explicatives sur les variables à expliquer. Les rayons des zones tampon dépendent de la résolution spatiale des données sources. Ils vont de 20 mètres (pour les métriques paysagères calculées avec l’images SPOT6 à 1.5 m de résolution spatiale) à 2 km (distance maximum de vol dans la vie d’un moustique selon la bibliographie). Voir le tableau ci-dessous, colonne “buffer_size_extraction_m”, pour le détail des tailles de buffer utilisées pour chaque variable.
Les variables temporelles ont été extraites durant les 120 jours (4 mois) précédent la date de capture ou d’enquête épidémiologique, afin de mesurer l’éventuelle influence temporelle de ces variables explicatives sur les variables à expliquer (cela sera fait via la génération de cross correlation maps, cf. section plan d’analyse de données). Cette durée maximale de 120 jours a été utilisée dans d’autres études similaires. Elle permet de couvrir largement la durée de vie d’un vecteur (en général, maximum 1 mois) et de prendre en compte l’éventuelle influence des générations de vecteurs précédent. Notons que les données temporelles ont des résolutions temporelles (i.e. temps de revisite) variables selon les sources : cela va du quotidien (températures de MODIS, précipitations de GPM) au mensuel (lumières noctures) en passant par l’hebdomadaire (évapotranspiration ou indices de végétation de MODIS), la dizaines de jours (Sentinel-1), etc.
Ainsi, si n est le nombre de jours précédent la capture pour lesquels nous avons extrait les données temporelles (ici n = 120), r la résolution temporelle de la donnée source de la variable à extraire, et m le nombre de zones tampon pour une variable spatiale considérée :
Chacune des variables issues de ces trois groupes fera l’objet d’un pré-traitement spécifique dans la modélisation (cf. section plan d’analyse de données).
code | short_name | long_name | type_group1 | type_group2 | type_group3 | spatial_variable | time_dependent_variable | max_time_lag | temporal_aggregation_days | unit | buffer_size_extraction_m | date_source_code | data_source_long_name | data_source_doi | spatial_resolution_m | revisit_time | latency_days | data_prepation_lineage | comments | table_db |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
TMIN1 | Daily minimum temperature | Daily minimum land surface temperature | Temperature | Medium temporal dynamicity | emergence | Yes | Yes | 120 | 1 | celcius degrees | 500, 1000, 2000 | MOD11A1.006, MYD11A1.006 | MODIS/Terra Land Surface Temperature/Emissivity Daily L3 Global 1km SIN Grid V006, MODIS/Aqua Land Surface Temperature/Emissivity Daily L3 Global 1km SIN Grid V006 | https://dx.doi.org/10.5067/MODIS/MOD11A1.006 https://dx.doi.org/10.5067/MODIS/MYD11A1.006 | 1000 | 1 | 2 |
|
NA | env_timeseries |
TMAX1 | Daily maximum temperature | Daily maximum land surface temperature | Temperature | Medium temporal dynamicity | emergence | Yes | Yes | 120 | 1 | celcius degrees | 500, 1000, 2000 | MOD11A1.006, MYD11A1.006 | MODIS/Terra Land Surface Temperature/Emissivity Daily L3 Global 1km SIN Grid V006, MODIS/Aqua Land Surface Temperature/Emissivity Daily L3 Global 1km SIN Grid V006 | https://dx.doi.org/10.5067/MODIS/MOD11A1.006 https://dx.doi.org/10.5067/MODIS/MYD11A1.006 | 1000 | 1 | 2 | Same as for TMIN1 but with the LST_Day_1km band and keeping the highest value available value for each pixel. | NA | env_timeseries |
TAMP1 | Daily temperature amplitude | Daily temperature amplitude | Temperature | Medium temporal dynamicity | emergence | Yes | Yes | 120 | 1 | celcius degrees | 500, 1000, 2000 | MOD11A1.006, MYD11A1.006 | MODIS/Terra Land Surface Temperature/Emissivity Daily L3 Global 1km SIN Grid V006, MODIS/Aqua Land Surface Temperature/Emissivity Daily L3 Global 1km SIN Grid V006 | https://dx.doi.org/10.5067/MODIS/MOD11A1.006 https://dx.doi.org/10.5067/MODIS/MYD11A1.006 | 1000 | 1 | 2 | Compute the difference between final TMAX1 and TMIN1 datasets. | NA | env_timeseries |
VNV8 | NDVI - 250 m | Average 8-days Normalized Difference Vegetation Index (250 m) | Land cover | Medium temporal dynamicity | emergence+dispersion | Yes | Yes | 120 | 8 | unitless | 500, 1000, 2000 | MOD13Q1.006, MYD13Q1.006 | MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid V006 MODIS/Aqua Vegetation Indices 16-Day L3 Global 250m SIN Grid V006 | https://dx.doi.org/10.5067/MODIS/MOD13Q1.006 https://dx.doi.org/10.5067/MODIS/MYD13Q1.006 | 250 | 8 | 20 |
|
The MODIS Vegetation products are produced over a 16 days time frame, hence there are very few missing values in the raw products. In addition, Terra and Aqua Vegetation products are phased, enabling to have a 8-days temporal resolution. | env_timeseries |
VEV8 | EVI - 250 m | Average 8-days Enhanced Vegetation Index (250 m) | Land cover | Medium temporal dynamicity | emergence+dispersion | Yes | Yes | 120 | 8 | unitless | 500, 1000, 2000 | MOD13Q1.006, MYD13Q1.006 | MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid V006 MODIS/Aqua Vegetation Indices 16-Day L3 Global 250m SIN Grid V006 | https://dx.doi.org/10.5067/MODIS/MOD13Q1.006 https://dx.doi.org/10.5067/MODIS/MYD13Q1.006 | 250 | 8 | 20 | Same as for VNV8 but with the _250m_16_days_EVI band. | The MODIS Vegetation products are produced over a 16 days time frame, hence there are very few missing values in the raw products. In addition, Terra and Aqua Vegetation products are phased, enabling to have a 8-days temporal resolution. | env_timeseries |
EVT8 | Evapotranspiration | Average 8-days evapotranspiration | Moisture proxy | Medium temporal dynamicity | emergence | Yes | Yes | 120 | 8 | kg/m�/8day | 500, 1000, 2000 | MOD16A2.006, MYD16A2.006 | MODIS/Terra Net Evapotranspiration 8-Day L4 Global 500m SIN Grid V006 MODIS/Aqua Net Evapotranspiration 8-Day L4 Global 500m SIN Grid V006 | https://dx.doi.org/10.5067/MODIS/MOD16A2.006 https://dx.doi.org/10.5067/MODIS/MYD16A2.006 | 500 | 8 | 20 |
|
The MODIS Evapotranspiration products are produced over a 8 days time frame, but are not phased (unlike MODIS Vegetation). This is the reason why we combine Terra and Aqua products to produce the final dataset. | env_timeseries |
SMO1 | Soil moisture | Daily soil moisture | Moisture proxy | Medium temporal dynamicity | emergence | Yes | Yes | 120 | 1 | mm3/mm3 | 500, 1000, 2000 | SPL3SMP_E.003 | SMAP Enhanced L3 Radiometer Global Daily 9 km EASE-Grid Soil Moisture, Version 3 | https://doi.org/10.5067/T90W6VRLCBHI | 9000 | 2/3 | 1 |
|
NA | env_timeseries |
RFD1_F | Rainfall | Daily total precipitation (research quality data) | Moisture proxy | Medium temporal dynamicity | emergence | Yes | Yes | 120 | 1 | mm | 500, 1000, 2000 | GPM_3IMERGDF.06 | GPM IMERG Final Precipitation L3 1 day 0.1 degree x 0.1 degree V06 | https://doi.org/10.5067/GPM/IMERGDF/DAY/06 | 10000 resampled to 250 | 1 | 1 | NA | NA | env_timeseries |
VNV30 | NDVI - 10 m | Average monthly Normalized Difference Vegetation Index (10 m) | Land cover | Medium temporal dynamicity | emergence+dispersion | Yes | Yes | 120 | 30 | unitless | 100, 250, 500, 500, 1000, 2000 | Sentinel-2 | Copernicus Sentinel-2 Data | Contains modified Copernicus Sentinel data 2020 for Sentinel data | 10 | 5 | 0 |
|
NA | env_timeseries |
VMV30 | MNDVI - 10 m | Average monthly Modified Normalized Difference Vegetation Index | Land cover | Medium temporal dynamicity | emergence+dispersion | Yes | Yes | 120 | 30 | unitless | 100, 250, 500, 500, 1000, 2000 | Sentinel-2 | Copernicus Sentinel-2 Data | Contains modified Copernicus Sentinel data 2020 for Sentinel data | 10 | 5 | 0 | Same as for VNV30 | NA | env_timeseries |
WNW30 | NDWI - 10 m | Average monthly Normalized Difference Water Index | Land cover | Medium temporal dynamicity | emergence | Yes | Yes | 120 | 30 | unitless | 100, 250, 500, 500, 1000, 2000 | Sentinel-2 | Copernicus Sentinel-2 Data | Contains modified Copernicus Sentinel data 2020 for Sentinel data | 10 | 5 | 0 | Same as for VNV30 | NA | env_timeseries |
BRI30 | BRI - 10 m | Average monthly Brillance Index | Land cover | Medium temporal dynamicity | dispersion+resistence | Yes | Yes | 120 | 30 | unitless | 100, 250, 500, 500, 1000, 2000 | Sentinel-2 | Copernicus Sentinel-2 Data | Contains modified Copernicus Sentinel data 2020 for Sentinel data | 10 | 5 | 0 | Same as for VNV30 | NA | env_timeseries |
WVV10 | VV - 10 m | Average 10-days VV band | Land cover | Medium temporal dynamicity | emergence | Yes | Yes | 120 | 10 | NA | 100, 250, 500, 500, 1000, 2000 | Sentinel-1 | Copernicus Sentinel-1 Data | Contains modified Copernicus Sentinel data 2020 for Sentinel data | 10 | 10 | 0 | NA | NA | env_timeseries |
WVH10 | VH - 10 m | Average 10-days VH band | Land cover | Medium temporal dynamicity | emergence | Yes | Yes | 120 | 10 | NA | 100, 250, 500, 500, 1000, 2000 | Sentinel-1 | Copernicus Sentinel-1 Data | Contains modified Copernicus Sentinel data 2020 for Sentinel data | 10 | 10 | 0 | NA | NA | env_timeseries |
LIG30 | Nightly radiance | Average monthly nightly radiance | Lights | Medium temporal dynamicity | activity | Yes | Yes | 120 | 30 | NanoWatt/cm2/sr | 500, 1000, 2000 | VIIRS_DNB_MONTH | Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band (DNB) | https://doi.org/10.5067/VIIRS/VNP46A1.001 | 30 | 30 | 30 |
|
NA | env_timeseries |
DTL7 | Daytime length | Weekly daytime length | Daytime length | Medium temporal dynamicity | emergence | No | Yes | 120 | 7 | hours | NA | geosphere | geosphere R package | https://CRAN.R-project.org/package=geosphere | NA | 1 | NA | Get the daylength at the sampling point location using the R package ‘geosphere’ | NA | env_timeseries |
RFH | Rainfall during the night of catch | Proportion of half-hours with positive precipitation for the whole duration of the HLC | Rainfall | Short temporal dynamicity | activity | No | No | NA | NA | % | NA | GPM_3IMERGDHH.06 | GPM IMERG Final Precipitation L3 Half Hourly 0.1 degree x 0.1 degree V06 | https://doi.org/10.5067/GPM/IMERG/3B-HH/06 | 10000 resampled to 250 | 0.02 | NA |
|
NA | env_nightcatch |
WDR | Wind direction during the night of catch | Mean wind direction during the night of catch | Wind | Short temporal dynamicity | activity | No | No | NA | NA | degrees (0 to 360) | NA | ERA5 | ERA5 hourly data on single levels from 1979 to present | https://doi.org/10.24381/cds.adbb2d47 | 27000 | 0.04 | NA |
|
NA | env_nightcatch |
WSP | Wind speed during the night of catch | Mean wind speed direction during the night of catch | Wind | Short temporal dynamicity | activity | No | No | NA | NA | m/s | NA | ERA5 | ERA5 hourly data on single levels from 1979 to present | https://doi.org/10.24381/cds.adbb2d47 | 27000 resampled to 250 | 0.04 | NA | See WDR | NA | env_nightcatch |
NMA | Atmoshperic pressure during the night of catch | Mean atmospheric pressure during the night of catch | Atmospheric pressure | Short temporal dynamicity | activity | No | No | NA | NA | NA | NA | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_nightcatch |
NMT | Temperature during the night of catch | Mean temperature during the night of catch | Temperature | Short temporal dynamicity | activity | No | No | NA | NA | NA | NA | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_nightcatch_postedecapture |
NML | Luminosity during the night of catch | Mean luminosity during the night of catch | Lights | Short temporal dynamicity | activity | No | No | NA | NA | NA | NA | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_nightcatch_postedecapture |
NMH | Humidity during the night of catch | Mean humidity during the night of catch | Humidity | Short temporal dynamicity | activity | No | No | NA | NA | NA | NA | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_nightcatch_postedecapture |
NDP | Dew point during the night of catch | Mean dew point during the night of catch | Dew point | Short temporal dynamicity | activity | No | No | NA | NA | NA | NA | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_nightcatch_postedecapture |
LMN | Apparent magnitude of the Moon | Apparent magnitude of the Moon | Lights | Short temporal dynamicity | activity | No | No | NA | NA | unitless | NA | MIRIADE | The Virtual Observatory Solar System Object Ephemeris Generator | NA | NA | 1 | NA |
|
NA | env_nightcatch |
TEL | Elevation | Average elevation | Topography | Long or nonexistent temporal dynamicity | NA | Yes | No | NA | NA | meters above the see | 500, 1000, 2000 | SRTMGL1.003 | Digital Elevation Model from the NASA Shuttle Radar Topography Mission Global 1 arc second | https://dx.doi.org/10.5067/MEASURES/SRTM/SRTMGL1.003 | 30 | NA | NA |
|
NA | env_staticbuffer |
TSL | Slope | Average slope | Topography | Long or nonexistent temporal dynamicity | emergence | Yes | No | NA | NA | % | 500, 1000, 2000 | SRTMGL1.003 | Digital Elevation Model from the NASA Shuttle Radar Topography Mission Global 1 arc second | https://dx.doi.org/10.5067/MEASURES/SRTM/SRTMGL1.003 | 30 | NA | NA | See TEL | NA | env_staticbuffer |
TAS | Aspect | Average aspect | Topography | Long or nonexistent temporal dynamicity | emergence | Yes | No | NA | NA | NA | 500, 1000, 2000 | SRTMGL1.003 | Digital Elevation Model from the NASA Shuttle Radar Topography Mission Global 1 arc second | https://dx.doi.org/10.5067/MEASURES/SRTM/SRTMGL1.003 | 30 | NA | NA | See TEL | NA | env_staticbuffer |
TCI | Terrain Classification Index | Average Terrain Classification Index | Topography | Long or nonexistent temporal dynamicity | emergence | Yes | No | NA | NA | unitless | 500, 1000, 2000 | SRTMGL1.003 | Digital Elevation Model from the NASA Shuttle Radar Topography Mission Global 1 arc second | https://dx.doi.org/10.5067/MEASURES/SRTM/SRTMGL1.003 | 30 | NA | NA | See TEL | NA | env_staticbuffer |
TWI | Topographic Wetness Index | Average Topographic Wetness Index | Topography | Long or nonexistent temporal dynamicity | emergence | Yes | No | NA | NA | unitless | 500, 1000, 2000 | SRTMGL1.003 | Digital Elevation Model from the NASA Shuttle Radar Topography Mission Global 1 arc second | https://dx.doi.org/10.5067/MEASURES/SRTM/SRTMGL1.003 | 30 | NA | NA | See TEL | NA | env_staticbuffer |
WAC | Accumulation | Average accumulation | Moisture proxy | Long or nonexistent temporal dynamicity | emergence | Yes | No | NA | NA | ha | 500, 1000, 2000 | SRTMGL1.003 | Digital Elevation Model from the NASA Shuttle Radar Topography Mission Global 1 arc second | https://dx.doi.org/10.5067/MEASURES/SRTM/SRTMGL1.003 | 30 | NA | NA | See TEL | NA | env_staticbuffer |
WAD | Distance to hydrographic network | Average distance to hydrographic network | Moisture proxy | Long or nonexistent temporal dynamicity | emergence | Yes | No | NA | NA | m | 500, 1000, 2000 | SRTMGL1.003 | Digital Elevation Model from the NASA Shuttle Radar Topography Mission Global 1 arc second | https://dx.doi.org/10.5067/MEASURES/SRTM/SRTMGL1.003 | NA | NA | NA |
|
NA | env_staticbuffer |
WMD | Shortest distance to hydrographic network | Shortest distance to theoretical hydrographic network | Moisture proxy | Long or nonexistent temporal dynamicity | emergence | No | No | NA | NA | m | 500, 1000, 2000 | SRTMGL1.003 | Digital Elevation Model from the NASA Shuttle Radar Topography Mission Global 1 arc second | https://dx.doi.org/10.5067/MEASURES/SRTM/SRTMGL1.003 | NA | NA | NA | See WAD | NA | env_staticnobuffer |
WLS | Length of the hydrographic network | Total length of the hydrographic network | Moisture proxy | Long or nonexistent temporal dynamicity | emergence | Yes | No | NA | NA | m | 500, 1000, 2000 | SRTMGL1.003 | Digital Elevation Model from the NASA Shuttle Radar Topography Mission Global 1 arc second | https://dx.doi.org/10.5067/MEASURES/SRTM/SRTMGL1.003 | NA | NA | NA | See WAD | NA | env_staticbuffer |
WAL | Accumulation / distance to sampling point | Ratio between the accumulation and the distance to sampling point | Moisture proxy | Long or nonexistent temporal dynamicity | emergence | Yes | No | NA | NA | ha/m | 500, 1000, 2000 | SRTMGL1.003 | Digital Elevation Model from the NASA Shuttle Radar Topography Mission Global 1 arc second | https://dx.doi.org/10.5067/MEASURES/SRTM/SRTMGL1.003 | NA | NA | NA | See WAD | NA | env_staticbuffer |
POP | Human population | Total human population (census/ground data) | Hosts | Long or nonexistent temporal dynamicity | reproduction | Yes | No | NA | NA | person | 100, 500 | REACT database | database of the REACT project | NA | household | NA | NA |
|
NA | env_staticbuffer |
BDE | Distance to the edge of the village | Distance to the nearest edge of the village | Land cover | Long or nonexistent temporal dynamicity | dispersion | No | No | NA | NA | m | NA | REACT database | database of the REACT project | NA | NA | NA | NA |
|
NA | env_staticnobuffer |
BCH | Degree of clustering of the households in the village | Degree of clustering or ordering of the households | Land cover | Long or nonexistent temporal dynamicity | dispersion | No | No | NA | NA | NA | NA | REACT database | database of the REACT project | NA | NA | NA | NA |
|
NA | env_staticbuffer |
HYS | Proportion of hydromorphic soils | Proportion of hydromorphic soils | Moisture proxy | Long or nonexistent temporal dynamicity | emergence | Yes | No | NA | NA | % | 100, 250, 500, 500, 1000, 2000 | IRD | Ressources en sols : carte à 1:500 000 des unités agronomiques déduites de la carte pédologique : ouest sud. (IN) Notice des cartes de ressources en sols de la Haute Volta Korhogo : carte des segments pédologiques. (IN) Etude pédologique de la région de Boundiali - Korhogo (Côte-d’Ivoire) : cartographie et typologie sommaire des sols : feuille Boundiali, feuille Korhogo à 1/200.000 | http://sphaera.cartographie.ird.fr/carte.php?num=2257&pays=BURKINA+FASO&iso=BFA http://sphaera.cartographie.ird.fr/carte.php?num=942&pays=COTE+D%2527IVOIRE&iso=CIV | NA | NA | NA |
|
NA | env_staticbuffer |
VCP | Vector control measure | Implementation of complementary vector control measures (other than LLIN) | Vector control | Vector control measures | longevity+resistence | No | No | NA | NA | 1 = presence ; 0 = absence | NA | REACT database | database of the REACT project | NA | NA | NA | NA |
|
NA | env_staticnobuffer |
VCM | Vector control measure type | Vector control measure type (LLIN / IRS / anti-larval fighting / communication strategies) | Vector control | Vector control measures | longevity+resistence | No | No | NA | NA | NA | NA | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_staticnobuffer |
VCT | Time since introduction of vector control measures | Numer of days since introduction of vector control measures (to study rebounding residual transmission, cf kileen) | Vector control | Vector control measures | longevity+resistence | No | No | NA | NA | number of days after the start of the intervention | NA | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_staticnobuffer |
IEH | Household interior / exterior | Interior or exterior of the household | Exophagy | Long or nonexistent temporal dynamicity | dispersion+resistence | No | No | NA | NA | NA | NA | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | NA |
LED | Land cover - edge density | Cumulated length of each land cover class perimeter (class level) | Land cover | Long or nonexistent temporal dynamicity | emergence+dispersion+activity | Yes | No | NA | NA | various units depending on the metric | 20, 50, 100, 250, 500, 500, 1000, 2000 | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_lsm |
LNP | Land cover - number of patches | Number of patches for each land cover class (class level) | Land cover | Long or nonexistent temporal dynamicity | emergence+dispersion+activity | Yes | No | NA | NA | various units depending on the metric | 20, 50, 100, 250, 500, 500, 1000, 2000 | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_lsm |
LPL | Land cover - percentage of landscape | Proportion of the total surface occupied by each land cover class (class level) | Land cover | Long or nonexistent temporal dynamicity | emergence+dispersion+activity | Yes | No | NA | NA | various units depending on the metric | 20, 50, 100, 250, 500, 500, 1000, 2000 | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_lsm |
LRD | Land cover - patch richness density | Number of different land cover classes (landscape level) | Land cover | Long or nonexistent temporal dynamicity | emergence+dispersion+activity | Yes | No | NA | NA | various units depending on the metric | 20, 50, 100, 250, 500, 500, 1000, 2000 | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_lsm |
LSD | Land cover - shannon’s diversity index | Shannon’s diversity (landscape level) | Land cover | Long or nonexistent temporal dynamicity | emergence+dispersion+activity | Yes | No | NA | NA | various units depending on the metric | 20, 50, 100, 250, 500, 500, 1000, 2000 | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_lsm |
LSE | Land cover - simspon’s evenness index | Simspon’s diversity (landscape level) | Land cover | Long or nonexistent temporal dynamicity | emergence+dispersion+activity | Yes | No | NA | NA | various units depending on the metric | 20, 50, 100, 250, 500, 500, 1000, 2000 | REACT database | database of the REACT project | NA | NA | NA | NA | NA | NA | env_lsm |
Les données que nous modélisons présentent les particularités suivantes :
Nous devrons prendre en compte ces particularités dans le choix des modèles à utiliser. L’analyse des résidus des modèles générés (incluant notamment l’analyse de l’autocorrelation spatiale) permettra de valider les modèles.
Nous présentons ici très succintement les avantages et limites de quelques modèles mathématiques usuels au regard des particularités présentées dans la section précédente :
Comme décrit précédemment, la variable à modéliser suit une distribution binomiale négative avec excès de zéros. Ces distributions sont souvent analysées avec des modèles de “haie”. Dans un modèle de haie, deux processus différents sont modélisés : les comptes nuls vs. les comptes positifs (modèle de présence/absence) et l’abondance des comptes positifs uniquement (modèle d’abondance). Cette séparation est motivée à la fois par des raisons statistiques (un seul et unique modèle gérerait mal ce type de distribution extrème) et par des considérations pratiques : les processus (biologiques, écologiques, etc.) conduisant à la présence des vecteurs peuvent être différents de ceux conduisant à leur abondance. Nous proposons ici de procéder de cette manière, en modélisant séparement la présence et l’abondance des vecteurs.
Nous établirons donc deux modèles : un modèle de présence/absence (classification) et un modèle d’abondance (regression). Comme l’on pose l’hypothèse que les processus conduisant à ces deux situations (présence/absence et abondance) sont différents, idéalement, chacun de ces deux modèles sera établi de manière indépendante, notamment dans la phase de séléction supervisée des variables explicatives à intégrer dans le modèle multivarié.
Les variables explicatives, en particulier les variables spatio-temporelles et les variables d’occupation du sol, sont très nombreuses et ainsi, quel que soit le modèle multivarié utilisé, il est nécéssaire de les filtrer avant leur intégration dans ce denier. Cette section présente les méthodes envisagées pour séléctionner les variables à insérer dans les modèles multivariés.
Comme décrit précédemment, nous avons collecté les variables spatio-temporelles sur une durée de 120 jours avant la date de capture. Plutôt que d’étudier la relation entre la variable à expliquer et la variable explicative pour chaque date prise séparemment, il est intéressant d’étudier des combinaisons d’intervalles temporels. Par exemple : il est plus probable que la somme des précipitations entre n jours et m jours avant la capture explique mieux la densité aggressive que les précipitations ponctuelles j jours avant la capture.
Pour cela, nous utilisons des cross-correlation maps (CCM) pour étudier les variables temporelles. Pour un lag de temps minimal et un lag maximal fixés ainsi qu’une fonction imposée (par exemple la moyenne ou la somme) les CCM déterminent les deux lags entre lesquels la fonction a la plus forte corrélation avec la variable étudiée. La valeur de la variable explicative entre ces deux lags temporels est ensuite retenue pour les modèles multivariés.
L’objectif ici est donc de définir une méthode pour corréler la variable à expliquer et chaque lag temporel des variables explicatives. Dans la littérature, ce sont des coefficients de corrélation classiques (pearson ou spearman) qui sont utilisés. Ces méthodes présentent l’avantage d’être rapides (élément non négligeable car les combinaisons de lags sont très nombreuses - jusqu’à 7500 pour les variables à résolution temporelle quotidienne) et de fournir à la fois le sens et la force de la corrélation, mais présentent les inconvénients suivants :
A titre d’exemple, lLa figure suivante donne la CCM pour la série “abondance tous Anopheles confondus” et la variable “précipitations quotidiennes” en utilisant un simple coefficient de Spearman :
Les questions auxquelles nous devons répondre quand à la séléction des variables explicatives temporelles sont donc les suivantes :
- Quel(s) test(s) statistique(s) implémenter pour séléctionner la(es) variable(s) à retenir dans le modèle multivarié de présence / absence ?
- Quel(s) test statistique(s) implémenter pour séléctionner la(es) variable(s) à retenir dans le modèle multivarié d’abondance ?
Pas encore d’idée ici…
Proposition n°1 : Voir code R, ligne 137)
Le package ‘correlation’ de R permet de calculer des coefficients de corrélations classiques (pearson, spearman, etc) tout en intégrant un facteur comme effet aléatoire (voir https://easystats.github.io/correlation/articles/multilevel.html ). Aussi, nous pourrions utiliser ce package pour calculer un coefficient de Spearman en intégrant le village en effet aléatoire. avantages : prise en compte de l’effet aléatoire au niveau du village, force et signe de l’effet, détéction de relation non-linéaires
inconvénients : non détéction de relations non-monotones
Proposition n°2 : Voir code R, ligne 146)
Avec le même package, on peut utiliser le coefficient de corrélation de distance. Ce coefficient permet non seulement de détécter les relation non-linéaires mais aussi les relations non-monotones (cf https://en.wikipedia.org/wiki/Distance_correlation). Par contre, il n’est pas possible d’intégrer un facteur comme effet aléatoire.
avantages : détéction de relation non-linéaires + non-monotones
inconvénients : non prise en compte de l’effet aléatoire village
Proposition n°3 : Voir code R, ligne 155)
Le package ‘glmertree’ (Generalized linear mixed-effects model trees) de R permet de générer des arbres de décision tout en intégrant des effets aléatoires comme dans un modèle linéaire mixte. “The estimation algorithm iterates between (1) estimation of the tree given an offset of random effects, and (2) estimation of the random effects given the tree structure.” L’intéret d’utiliser un arbre de décision est sa capacité à détecter des relations complexes entre les variables. Par ailleurs l’introduction d’effets aléatoires permet de gérer la problématique de non indépendance des données. Finalement, un des avantages des glmertree est la gestion de tous types de variables réponses (quantitatives continues, quantitatives discrètes, binomiales).
Voir help(glmertree) section Details, pour un court résumé de la théorie des Generalized linear mixed-effects model trees ; https://link.springer.com/article/10.3758/s13428-017-0971-x pour l’article scientifique complet ; et vignette(“glmertree”) pour un tutoriel.
Cependant la question dans cette méthode est celle de la valeur quantitative finale à considérer dans les CCM. Dans les propositions 1 et 2 (voir ci-dessus), c’est le coefficient de corrélation qui est représenté, et la série présentant le coefficient le plus important est séléctionnée pour la suite des travaux. Dans le cas des glmertree, quelle valeur considérer ?
Voir code R, ligne 195)
Les variables d’occupation du sol sont des métriques paysagères (landscape metrics) qui ont été calculées selon plusieurs rayons de zones tampon autour de chaque point de capture (pour prendre en compte l’influence de la taille de la zone tampon considérée).
Ainsi, nous avons au total environ 300 variables d’occupation du sol, qu’il est nécéssaire de filtrer avant intégration dans les modèles multivariés. Les problématiques sont les mêmes que pour les variables temporelles :
- comment séléctionner les variables à intégrer dans les modèles multivariés de présence/absence ?
- comment séléctionner les variables à intégrer dans les modèles multivariés d’abondance ?
à définir
à définir
Comme la proposition n°1 des variables spatio-temporelles, mais en prenant le numéro de mission comme facteur aléatoire.
Voir code R, ligne 209 à 244)
Les autres variables sont principalement des variables liées à la topographie (altitude, pente, réseau hydrographique théorique, etc.), à la démographie (population), aux conditions micro-climatiques la nuit de capture (vent, précipitations, température, humidité, etc.).
Si l’on utilise des modèles algorithmiques (type random forest) pour la modélisation multivariée, il ne sera pas nécéssaire d’effectuer une pré-séléction de ces variables (ces modèles gèrent efficacement un grand nombre de variables explicatives). Si l’on utilise un modèle paramétrique plus classique (GLMM, GAMM), on devra les filter.
Quel que soit le modèle multivarié utilisé, il est important de prendre en compte la multicollinéarité entre les variables explicatives. La présence de variables collinéaires peut impacter la qualité du modèle lui même et / ou l’interprétation (par exemple dans le cas des variables d’importance ou des partial dependence plots).
Idéalement l’identification de la multicollinéarité dépend du modèle utilisé en modélisation multivarié :
from https://explained.ai/rf-importance/
“A way to identify if a feature, x, is dependent on other features is to train a model using x as a dependent variable and all other features as independent variables (this is called Multicollinearity). Because random forests give us an easy out of bag error estimate, the feature dependence functions in rfpimp rely on random forest models. The R^2 prediction error from the model indicates how easy it is to predict feature x using the other features. The higher the score, the more dependent feature x is. The feature importance of non-x features predicting x give an indication of which features have predictive power for feature x”
Comme décrit précédemment, il semble intéressant dans notre cas d’utiliser des modèles algorithmiques (non paramétriques), de part la quantité importante de variables explicatives et la forte probabilité de présence de relations non-linéaires et intéractions complexes entre variables. Ces modèles, réputés “black box”, sont aujourd’hui interprétables gràce aux divers outils développés (voir https://christophm.github.io/interpretable-ml-book/agnostic.html).
Nous proposons ici de tenter d’utiliser les modèles de foret aléatoire (random forest) à la fois pour la modélisation de la présence/absence de vecteurs (modèle de classification) et de l’abondance (modèle de régression).
Cependant les modèles algorithmiques présentent des limites qu’il est important de considérer :
Ces limites ne sont pas un problème dans les méthodes statistiques “classiques” :
Comme nous allons le voir ces limites ont été identifiées et étudiées ces dernières années dans le cas des modèles algorithmiques et des solutions ont été proposées pour y remédier.
Solution n°1 : intégrer dans le modèle des variables renseigant la dimension géographique des données. Certaines études montrent que l’insertion de variables géographiques permettent de considérer l’auto-corrélation spatiale. La théorie est la suivante : si le phénomène étudié exhibe un caractère spatial qui a une importance, les variables géographiques seront prises en compte dans le modèle. Certaines études préconisent d’intégrer les coordonnées géographique des points observations (X et Y), d’autres la matrice de distance entre tous les points d’observation (https://peerj.com/articles/5518/), d’autres encore les coordonnées géographiques ainsi que les distances à certains points fixes de la zone d’étude (par exemple les 4 extrémités et le centre) (https://doi.org/10.1111%2Fejss.12687).
Solution n°2 : insérer des variables renseignant les effets aléatoires. Certaines études montrent que la simple introduction des effets aléatoires en tant que variables explicatives dans le modèle permet de considérer ces effets. Ici nous pourrions donc introduire les colonnes “codevillage” (e.g. ‘BOH’), “pointdecapture” (e.g. ‘1’) et “idpointdecapture” (e.g. ‘BOH1’) pour considérer les effets aléatoires.
Solution n°3 : mixed - effects random forest (MERF) : L’idée des mixed - effects random forest est d’insérer des effets aléatoires dans le modèle de random forest. Une explication des mixed-effect random forest est donnée ici : https://towardsdatascience.com/mixed-effects-random-forests-6ecbb85cb177 . Problème : les tentatives de constructions de packages de Mixed effect random forest sur R n’ont pas l’air très concluantes… Sur Python il existe un package qui semble assez aboutit (https://pypi.org/project/merf/) mais je ne maitrise pas Python.
Pour la classification, il s’agit de modéliser 2 classes : présence ou absence de piqure sur chaque nuit de capture. Pour toutes les espèce, la classe “absence” est largement majoritaire sur la classe “présence”, autrement dit, il y a eu beaucoup plus de nuits de capture avec absence de capture qu’avec présence d’au moins 1 capture. Aussi, nous sommes dans un cas typique d’“imbalanced data”, où les modèles algorithmiques performent mal.
Pour la régression, nous faisons face à un problèle similaire : bien que nous ayions enlevé la valeur 0 piqures grâce à la séparation de la modélisation “présence/absence” et “abondance”, la distribution de la série statistique reste binomiale négative : la valeur “1 piqure” concentre environ la moitié des valeurs, qui vont jusqu’à plus de quarante piqures parfois.
Des méthodes ont été proposées pour remédier à ce problème. En général, il s’agit de méthode de rééchantillonage des observations (augmentation du nombre d’observations de la classe minoritaire, ou diminution du nombre d’observations de la classe majoritaire, ou un équilibre entre les deux méthodes).
Solution 1 : rééchantilloner la(es) classe(s) minoritaire(s) (afin de ne pas perdre d’information concernant la classe majoritaire) avec la méthode SMOTE (Synthetic Minority Over-sampling Technique). En sortie de cet algorithme, le nombre d’observations par classe est équilibré.
Solution 2 : Une autre solution qui pourrait être testée consiste à pondérer les observations, et insérer ces poids dans le modèle. A priori, dans un modèle de foret aléatoire les observations à utiliser en entrée de chaque arbre de décision sont séléctionnées de manière aléatoire. Si l’on pondère les observations, une observation avec un poids important sera plus souvent séléctionnée pour générer les arbres de décision, qu’une observation avec un poids faible.
Solution 3 (classification binomiale uniquement) : les random forests quantile classifier (http://web.ccs.miami.edu/~hishwaran/papers/OI_PatternRecognition_2019.pdf) implémentés dans le package R rfsrc (fonction imbalanced.rfsrc)
[section à terminer]
Le modèle sera cross-validée en utilisant un méthode leave-one-mission-out : toutes les missions sauf une seront tour à tour utilisées pour entrainer le modèle, puis la mission qui n’a pas été utilisée pour l’entrainement sera utilisée pour la validation.
L’auto-corrélation des résidus du modèle (residual spatial autocorrelation) sera évaluée (via des variogrammes ou le calcul du I de Moran) pour le valider.
[section à terminer]
-> on mesure l’importance des variables (par ex. minimum depth et permutation importance)
-> on calcule les intéractions entre les variables les plus importantes
-> on mesure l’effet d’une variable explicative sur la variable réponse (avec des PDP et des ICE) (main effects) en prenant en compte l’ensemble des effets des autres variables explicatives et intéractions entre variables
-> on mesure l’intéraction entre 2 variables explicatives sur la variable réponse (avec des PDP) (second order effects), de même, en prenant en compte l’ensemble des effets des autres variables explicatives et intéractions entre variables
Comment interpréter les modèles ?
Dans le cas des modèles linéaires, les équations en sortie du modèle présentent l’intéret de procurer, pour chaque variable ou interactions entre variables : i) l’effet (sens et ampleur de l’influence grâce à resp. le signe et taille du coefficient), ii) la significativité de cet effet (grâce à la p-value). Les modèles algorithmiques complexes ne produisent pas d’équations comme celles-ci. Cependant, un certain nombre d’outils a été développé pour rendre ces modèles davantage intérprétables et ainsi extraire des informations similaires à celles que procurent les équations en sortie des modèles linéaires : les variables importantes, l’effet d’une variable ou d’une intéraction entre plusieurs variables sur la variable réponse (partial depence plots), etc. Ces méthodes sont pour la pluspart indépendantes du modèle qui est utilisé (model-agnostic interpretation methods).
Ainsi, si l’on utilise ces outils d’interprétation, les éléments les plus importants à considérer dans le choix du modèle à utiliser sont la capacité prédictive et évidemment le respect des hypothèses liées au modèle.