Modélisation de répartition

Bonjour,

Je suis étudiant en classe préparatoire et je réalise une étude sur la répartition des amphibiens et reptiles de Normandie. Nous cherchons à partir des observations et des données climatiques de Météo France à réaliser une carte des probabilités de présence. Ce travail a déjà été réaliser pour les reptiles de Suisse. A l'aide d'un logiciel SIG, nous avons pu construire un tableau des observations croisé avec les données climatiques (températures et pluviométries). Ainsi pour chaque espèce et pour maille de la Normandie nous avons le tableau qui indique par colonne les valeurs des variables climatiques et l'absence ou présence de l'espèce. Il nous reste à construire le modèle. Nous devons donc à partir du logiciel R avec la fonction modèle linéaire généralisé, de loi binomiale et de lien logit construire nos modèles. Nous nous étions orienté au début vers l'utilisation de tests statiques tel que le rapport de vraisemblance pour analyser nos modèles. Mais il est apparu que le critère AIC est plus simple d'utilisation et tout aussi efficace. Ce critère tient compte de l'exactitude du modèle et du nombre de variables à la fois. Mais nous ne savons comment procéder pour réaliser ces modèles, l'entrée des variables (et même leur ordre?) étant sensible. Quelqu'un saurait-il donc comment s'y prendre pour rentrer les variables dans le modèle, ou des documents l'expliquant? Nous savons qu'il faut tenir compte du principe de parcimonie et nous avons déjà établi nos hypothèses. Mais concrêtement, je n'ai pas trouvé si l'on devait rentrer d'abord toutes les variables puis éliminer petit à petit les variables les moins significatives pour voir si l'AIC diminue? Je sais qu'il faut vérifier l'hypothèse nulle, et qu'il faudra faire la même démarche sur un second jeu de données, mais je ne sais pas qu'elle est la méthode rigoureuse pour établir un modèle.
Nous devons faire le travail samedi. C'est une demande de dernière minute, mais s'il n'est pas trop tard merci d'avance!

Quentin

salut Quentin,

je vais essayer de t'aider mais je ne suis pas Statisticien, donc il se peut que je comprenne mal certains points.
En effet le critère d'Akaike (AIC) est le meilleur outil à utiliser dans ton cas selon moi. Il correspond à un Log de Max de Vraisemblance + 2ddl

Le test de modèle peut se faire en backward ou forward ; c'est a dire en partant du modele le plus complet que tu ais et auquel tu retires chacune des variables alternativement pour voir si tu peux gagner en vraissemblance, ou à l'inverse en partant du modele nul auquel tu tentes d'ajouter une des variables...
lis ce post qui explique cette idée:
https://ladoua.exprimetoi.net/biostatistiques-f59/forward-or-backward-t316.htm

selon moi tu devrais choisir le backward, si ton modele complet est bien ajusté à tes données ! Tu peux aussi faire les 2 (R fait tous, donc c'est pas lourd ) et voir si tu tombes sur la même chose ?...

Le plus simple sous R est de :
#rentrer ton modele complet
lmX=lm(Y~X1*X2*X3)
#test en backward, il partira du modele et cherchera le plus vraisemblable en retirant une a une les variables
Il n'y a pas de probleme d'ordre d'entrée, mais il testera d'abord tous les modeles en retirant alternativement une des interactions du plus haut niveau, puis de moindre niveau s'il y'a mieux, puis tes variables.
Le probleme est que si une interaction ne peut etre retirée, il s'arretera sans tester tes variables...
step(lmX)
# tu peux tester en forward, ou directement lui demander les deux
step(lmX,forward)
step(lmX,both)

tiens nous au courant de l'avancé de ton test,
en espérant t'avoir aidé...

J'y réfléchis! Merci beaucoup!
Je vous tiens au courant bien sûr!

Quentin

Bonjour,

Merci de votre réponse, qui amène à d’autres questions !
Tout d’abord pour l’entrée du modèle, nous prendrons toutes les variables sous formes de monômes d’ordre 1 de type glmX=glm(Y~X1*X2*X3) et non de pas avec des termes aux carrés, au cube…
Il faut ensuite choisir les variables d’entrée. Nous avons fait des hypothèses quant aux variables les plus vraisemblablement significatives. Rigoureusement faut-il rentrer toutes les variables (en éliminant les variables redondantes ou corrélées), ou directement les variables de notre présélection ? Dans ce dernier si nos hypothèses omettent une variable significative nous ne pourrons le savoir, mais dans l’autre cas nous savons que l’AIC a tendance à complexifier les modèles. Si l’on rentre la fonction step(glmX,forward) R nous sortira directement le modèle qu’il a trouvé le meilleur en utilisant le critère AIC ? Dans le cas de la méthode par backward, il nous sortirait donc juste l’AIC du dernier modèle retenu ? L’AIC étant un critère de comparaison entre modèles, avoir un seul AIC n’a donc pas de signification. On pourrait alors en entrant ensuite une autre formule avec d’autres variables (par exemple en remplaçant certaines variables redondantes pour voir laquelle est la plus significative), obtenir un autre modèle retenu avec son AIC et le comparer alors au précédent ? Et avec la fonction both, cela nous donne deux résultats, celui en backward et celui en forward? Je ne sais ce que signifie <none>, l’hypothèse nulle ?
On m’a dit aussi qu’il fallait absolument avant de commencer, vérifier l’hypothèse nulle. Mais concrètement comment le fait-on ? Nous sommes obligé de rentrer des paramètres dans le modèle, mais lesquels alors? L’intérêt serait donc de vérifier que l’AIC du modèle nul est bien supérieur à celui du modèle complet et du modèle retenu ? Et merci pour la fonction step, je ne l'avais encore jamais vue!

Désolé pour toutes ces questions !

Merci

Quentin

Salut,

Alors point par point...

Tout d’abord pour l’entrée du modèle, nous prendrons toutes les variables sous formes de monômes d’ordre 1 de type glmX=glm(Y~X1*X2*X3) et non de pas avec des termes aux carrés, au cube…

Dans un premier temps, je pense que c'est préférable, tu pourras complexifier si ton modele complet n'est pas correctement ajusté.

Il faut ensuite choisir les variables d’entrée. Nous avons fait des hypothèses quant aux variables les plus vraisemblablement significatives. Rigoureusement faut-il rentrer toutes les variables (en éliminant les variables redondantes ou corrélées), ou directement les variables de notre présélection ? Dans ce dernier si nos hypothèses omettent une variable significative nous ne pourrons le savoir, mais dans l’autre cas nous savons que l’AIC a tendance à complexifier les modèles.

En effet, il est tres important de ne pas laisser des variables corrélés dans ton modèle complet, choisi celle qui apporte le plus d'information (statistique) à ton modele. Ensuite, pour moi (de nombreux debats entre statisticiens existent sur ce sujet), il faut que tu entres toutes tes variables (relativement judicieuse tout de meme, tu ne mets pas le sexe de l'experimentateur pour tester la croissance d'une plante... quoi que ! ). Le probleme viendra peut être des interactions de plus haut niveau, s'il trouve une interaction d'ordre quatre significative, par exemple, il ne testera pas les interactions d'ordre 3 qui concerne ces 4 variables. Cela pourrait être ennuyeux... Nous verrons si c'est la cas, je ne sais pas trop comment agir spontannement face à ce cas à envisager (quelqu'un à une suggestion ??). Peut-etre te suggèrerais-je de te limiter , dans ce cas à des interactions d'ordre 2 !

Si l’on rentre la fonction step(glmX,forward) R nous sortira directement le modèle qu’il a trouvé le meilleur en utilisant le critère AIC ? Dans le cas de la méthode par backward, il nous sortirait donc juste l’AIC du dernier modèle retenu ?

Non, ce sont juste deux methodes différentes ayant le meme but trouver le meilleur modele !
Forward part de ton modele nul et essaie de rajouter alternativement chacune de tes variables. Il choisira le modele avec la variable donnant le plus faible AIC, puis testera si l'ajout d'une seconde peut la diminuer encore, ainsi de suite
Backward (que je preconise et qui est l'argument par defaut sous R)
procede a l'inverse. Il part de ton complet et essaie d'enlever chacun des variables pour voir si tu peux diminuer l'AIC

L’AIC étant un critère de comparaison entre modèles, avoir un seul AIC n’a donc pas de signification.

En effet !

On pourrait alors en entrant ensuite une autre formule avec d’autres variables (par exemple en remplaçant certaines variables redondantes pour voir laquelle est la plus significative), obtenir un autre modèle retenu avec son AIC et le comparer alors au précédent ?

euh... je ne pense pas que ce swoit très judicieux. Si tu choisis ta variable considéré comme la plus explicative, ce serait du bricolage que de la changer (tu peux le faire pour le coté informatif, de manière informelle, mais c'est tout, a mon sens)

Et avec la fonction both, cela nous donne deux résultats, celui en backward et celui en forward? Je ne sais ce que signifie , l’hypothèse nulle ?

Je n'ai jamais utilisé "both" mais c'est ce que je dirais spontannement.
Regarde l'aide dans R, pour mieux connaitre la fonction :
>?step
pour le , je pense que tu parle de la sortie de step, du genre :
Df Sum of Sq RSS AIC
- Agriculture:Examination 1 1.9 2836.6 204.7
2834.6 206.7
- Agriculture:Education 1 258.8 3093.5 208.8
- Examination:Education 1 340.6 3175.3 210.0

dans ce cas, <none> correspond à ton modele complet (ou retenu si tu es en deuxieme etape...) qui contient toute tes variables. Les autres correspondent à ce modèle -la variable (ex/ premiere ligne, mon modele moins l'intéraction Agric:Exam est meilleur que mon complet, il sera retenu et correspondra au <none> de l'etape suivant.)
Le test s'arrete quand le <none> est le meilleur modele, soit qu'on ne peut plus retirer (ou ajouter en forward) de variable.

On m’a dit aussi qu’il fallait absolument avant de commencer, vérifier l’hypothèse nulle. Mais concrètement comment le fait-on ?

Je ne connaissais pas cette précaution, et ne voit pas ce qu'elle pourrait signifier... en revanche il faut ABSOLUEMENT t'assurer que ton modele complet s'ajuste bien à tes données, qu'il les explique suffisamment. Sinon, cela signifie qu'il te manque quelques choses et que tu n'as pas le vrai modele complet !...
Pour le cas où le modele nul ait l'AIC la plus faible, c'est rien ne peut expliquer ton jeu de données (c'est rare) mais pour moi ce n'est pas un réel problème de méthode, donc, désolé, je ne vois pas de quoi tu parles:P

Bon courage

Adrien

PS : pour tes prochains messages, essaie de sauter des lignes pour gagner en visibilité

Bonjour,

Désolé pour la typographie !

Merci beaucoup, grâce à vos explications, nous avons pu manipuler step et mieux cerner ce que nous devions faire! Pour le modèle nul, c'est ce que l'on m'avait dit, mais je rejoignais plutôt votre idée. Il nous reste finalement quelques autres détails à régler pour nos matrices présence/absence et paramètres. Nous avons fait une sélection dans toutes les variables de météo France, selon les conseils avisés d’un excellent herpétologue.

Si cela vous intéresse, je vous présente rapidement notre démarche. Nous avons quelques données bibliographiques pour la répartition du lézard vert par exemple dont la limite correspondrait à l’isotherme de 18°C de Juillet et une pluviométrie annuelle inférieure à 600mm. En réalité ces données ne collent d’évidence pas pour la Normandie, et elles ont de plus été trouvée par des superpositions de carte et non des modélisation. C’est pourquoi notre but est de trouver quels sont les paramètres qui vont influer sur la répartition. Nous avons fait une sélection de paramètres, numéroté 01, …, 12 pour le mois (janvier, …, décembre) et 13 pour l’année.

La pluviométrie dans les études est considérée à l’année, on ne peut pas la considérer au mois étant donné ses variations déjà importantes d’une année à l’autre. Nous avons décider de la considérer par saison et par an. Nous prendrons les saisons au sens météorologique :
Printemps : mars, avril, mai
Eté : juin, juillet, août
Automne : septembre, octobre, novembre
Hiver : décembre, janvier, février

Voici les paramètres utilisés :

rrmoy01,…, rrmoy13 : hauteur de précipitations (en mm) par saison et annuelle
njrr01,…, njrr13 : nombre de jours avec hauteur de précipitation >= 1 mm par saison et annuel
tnmoy01,…, tnmoy13 : température minimale (en °C) mensuelle et annuelle
txmoy01,…, txmoy13 : température maximale (en °C) mensuelle et annuelle
njtx2513 : nombre de jours avec température maximale >= 25 °C annuel
njgel13 : nombre de jours avec gelée sous abri mensuel et annuel

Le nombre de paramètres que nous allons rentré est peut-être important ? mais nous ne savons pas lesquels se dégagerons exactement même si nous avons quelques idées dessus. Nous avons éliminer certains paramêtres qui étaient corrélés à ceux-ci. Peut-être y en a t'il encore qui sont corrélés, mais par exemple le nombre de jours avec température maximale >= 25 °C annuel ne fournit en fait pas la même information que la température maximale (en °C) annuelle (à la vue des cartes sous SIG), et des personnes discutent de la plus grande pertinence de l'un ou de l'autre. C'est pour cela que nous pensions garder des variables proches pour que le modèle nous sorte la plus pertinente. Peut-être y a-il là encore un problème de corrélation?

Pour les amphibiens-reptiles, les saisons météorologiques correspondent globalement à leur cycle biologique. Le printemps sera important pour la sortie d’hivernage et la reproduction, l’été pour ses éventuels périodes de stress hydrique, l’automne pour l’éclosion, et l’hiver ne devant pas être trop rude. De plus la pluviométrie est importante pour les amphibiens l’hiver puisque c’est pendant cette période que les mares se remplissent. C’est pour cela que tout a potentiellement une importance !

Nous pensons donc établir un premier modèle comme suit pour les amphibiens ( pour les reptiles nous retirerons la pluviométrie l’hiver) :

P(Y) = njgel_année + njrr_hiver+ njrr_printemps+ njrr_été+ njrr_automne + njrr_année + njtx25année + rrmoy_hiver + rrmoy_printemps + rrmoy_été + rrmoy_automne + rrmoy_année + tnmoyh_hiver + tnmoy_printemps + tnmoy_été + tnmoy_automne + txmoy_hiver + txmoy_printemps + txmoy_été + txmoy_automne + tnmoy_année + txmoy_année

Du premier modèle obtenu nous souhaiterions, aux couples paramètre-saison sélectionner ajouté les couples paramètre-mois de la saison correspondante. Ce second modèle sélectionnerait un mois de la saison si c’est celui-ci qui est significatif, ou conserverait le couple paramètre-saison si c’est à l’échelle de la saison que le paramètre prend tout son poids.
On pourrait obtenir une précision plus fine. Cette approche n’est sûrement pas très rigoureuse mais à notre niveau l’approximation est peut-être négligeable ? Parce que rentrer tous les couples paramètre-mois donnerait un nombre considérable de variables, or on sait que l’AIC a tendance à complexifier les modèles. L’échelle des saisons est très intéressante, mais nous souhaiterions aussi vérifier les valeurs bibliographiques qui elles sont au mois (peut-être sans grande signification d’ailleurs !). Voici donc la démarche à laquelle nous pensions, pas très rigoureuse mais essayant à notre niveau de répondre à nos préoccupations.

Nous allons aussi faire un modèle uniquement avec les indices d’aridités par saison et annuel que nous allons calculer avec les données que nous avons, ce qui nous permettra également de vérifier les valeurs bibliographiques dessus.

Je ne sais pas ce que vous en penser ?

Quentin

Salut,

Premiere chose, on peut se tutoyer je pense, nous sommes sur un forum d'étudiant... ça sera plus sympa
Deux questions par rapport à ce que tu nous propose.

1- pourquoi n'utilises tu pas une variable "mois" plutot que de créer 12 variables a chaque fois ? (l'annuel reviendra a la somme de ta variable), et tu ajoutes cette variable confondante à ton modele et cela diminue par 13 ton nb de variable...

2- quelle est ta question (hypothese), estceque toutes les variables que tu proposes sont potentiellement des explicatives d'interet ou y'a til des effets confondants que tu ajoutes pour ne pas perdre de l'info (ex: année de donnée...)

3- (c'est pas une question mais une remarque) Tu proposes un modele de type additif ! Je pense qu'au moins ajouter les interactions d'ordre2 entre tes variables d'interet est indispensable... il se peut qu'elle impacte les unes sur les autres.
exemple, l'impact du nombre de jour de gel diffère en fonction de la temp min obs, ce n'est donc pas simplement des effets additifs mais interactifs.

4- Effets corrélés conservés, c'est mal ! dans ce cas, je te préconise de faire deux modèles avec chacune de ces deux variables tour a tour, et de voir le plus vraisemblable. mais garder les deux risque de fortement biaisé le step.

Voila, affaire a suivre...

Bonjour,

il faut tout de même d'abord avoir une idée du nombre de sites considérés afin de savoir si le modèle peut correctement ajuster quelque chose. Un plan expérimental équilibré est une base à un "bon" modèle. Plus on augmente le nombre de variable explicative plus il faut de site notamment si on ajoute les interactions. En plus il faut se garder un jeu de validation car il n'est pas certains, surtout avec une procédure stepwise, que le modèle produit soit le meilleur modèle pour prédire l'occurence des espèces.
Personnellement, j'évite aussi d'utiliser des variables catégorielles (en classe) pour prédire une variable binaire. Je trouve que ces varaibles ont une facheuse tendance à trop bien ajuster des données elle aussi en classe (l'occurrence est en classe : oui ou non, 0 ou 1).

D'autre part l'AIC n'est pas le seul critère. Tu peux aussi utiliser le BIC qui est l'équivalent bayésien de l'AIC.

Pour être réellement rigoureux, je pense aussi qu'il vaut mieux établir tous les modèles possibles et leur AIC (ou BIC) associé pour voir quel est le meilleur modèle. Par ailleurs, même en utilisant l'AIC on n'a pas forcément le meilleur modèle, il faut toujours passer par un examen approfondi des résidus du modèle afin de savoir si le modèle ajuste correctement les données.

Bon courage.

Nicolas

Oui on peut se tutoyer en effet! C’est que je n’ai pas trop l’habitude de tutoyer à distance.

1. Je n’ai pas très bien compris pour la création d’une variable ‘mois’. Parce que la pluviométrie d’un mois peut être significative sans que ce soit le cas de ses températures min et max.

2. Mes hypothèses sont que si par exemple, la pluviométrie est très significative à l’année peut-être n’est-elle en réalité très significative que en été (variable confondante donc). Mais c’est vrai que c’est une bonne question, et peut-on y répondre directement ? soit que ce soit la pluviométrie à l’année et celle de l’été qui sont des variables explicatives (il faudrait un certain cumul de pluie à l’année pour nos espèces et un période de stress hydrique plus ou moins importante par ex) ou soit qu’elles sont des variables confondantes (dans ce cas la pluviométrie printemps-automne-hiver serait globalement une constante en Normandie ou tout du moins la variation n’aurait pas beaucoup d’influence sur nos espèces).
Donc en fait je partais au départ pour des variables saisons/année confondantes mais finalement après réflexion ta question me pose un problème…

3. Je ne sais pas du tout faire plus qu’un modèle additif et en particulier introduire des interaction d’ordre 2 !

4. J’ai fait une erreur en écrivant, je n’ai le nombre de jours avec gelées qu’à l’année. Donc pour l’hiver et le printemps Tmin et njgel13 serait plus ou moins liés mais ils ne donnent pas la même information. Il peut très bien n’y avoir qu’une journée au printemps de gelée… Donc je verrais plutôt ces variables comme deux variables explicatives.

Mais le step va essayer de retirer l’une des variables après l’autre, donc celle qui sera la moins pertinente sera quand même éliminée ?

Pour l’histoire de mes variables en saison et en mois, j’ai discuter aujourd’hui avec un autre herpétologue et il me confirme que les variables prennent sens pour des périodes et non pour des mois, et par conséquent pour nous simplifier la tâche on pourrait très bien se contenter des saisons sans aller jusqu’au mois?

En fait nous n’avons pas énormément d’observations, il faut que je regarde exactement mais nous en avons entre 100 et 200 selon les espèces. De plus je n’ai reçu aucune formation en statistique, ce que je sais sur les stats je l’ai appris cet année pour pouvoir mettre à bien ce travail. Nous souhaitons donc faire cette approche de modélisation, de manière à obtenir un résultat exploitable, sans faire de grossière erreur de conceptualisation, mais en sachant que nous ne pourrons exploiter au maximum des possibilités nos données. Nous avons donc conscience qu’un jeu de validation serait nécessaire mais avec le nombre de données que l’on a et notre qualification nous allons sûrement laisser de côté cette étape qui serait pourtant de rigueur.
De même je sais en effet qu’il existent d’autres critères, et sûrement mieux que l’AIC mais celui-ci déjà sous R et donc plus simple d’utilisation pour nous. Nous regarderons aussi sûrement l’exactitude des modèles qualitativement sous SIG en superposant nos cartes de probabilités de présence et de répartition réelle. Ces modèles doivent nous donner une idée de l’aire climatique potentiellement favorable à nos espèces. Car en plus nous savons que nos espèces dépendent d’autres facteurs que les facteurs climatiques, et certaines même ont une répartition probablement plus déterminée par le sol ou le paysage. Nous pourrions ainsi plus ou moins valider l’hypothèse selon laquelle ces dernières espèces dépendent davantage d’autres facteurs, en obtenant des aires potentielles qui ne se superposeraient pas du tout.

Merci de votre intérêt !

Quentin

Il n'est pas évident de s'étaler sur des explications aussi larges ici. Cela avancerait bien plus vite sur une discussion directe...

Est-ce-que ton projet justifierait le fait que tu viennes nous rendre visite sur Lyon pour discuter de tout ça autour d'un verre ?

C'est vrai que ce n'est pas facile d'en parler comme ça, et cela doit être assez frustrant pour toi vu mes connaissances...

Ta proposition est super sympa! mais malheureusement ce n'est pas possible pour moi parce que je suis en concours dans pas très longtemps...

Lorsque je fais quelquechose je le fais toujours à fond, et en plus j'ai vraiment accroché pour ce travail. Mais le niveau d'exigence demandé est sûrement bien moindre et je passe déjà beaucoup de temps dessus en négligeant d'autres matières. Déjà que je n'avais jamais vu de stats! Je vais essayé de faire du mieux que je peux, en ayant conscience des difficultés, mais je pense me limiter: comme il est avéré que les mois n'ont pas de signification biologique, je vais donc travailler avec les paramètres-saisons. Et puis je pourrais aussi faire un modèle avec les paramètres-saisons, un autre avec les paramètres-année, et un troisième avec l'ensemble. En comparant les AIC je pourrais ainsi voir quel sera le meilleur modèle, en prenant en compte la difficulté liée aux saisons et à l'année?

Quentin

Nous avons enfin fini nos manip et nos modélisations. Il y a pas mal de choses intéressantes. Finalement nous avons tenu compte des corrélations avec la fonction cor de R. Je vais rentrer en période de concours donc je n'ai pas trop le temps en ce moment, mais si tu veux je peux t'envoyer les résultats par mail?

Quentin

Je ne dis pas non !

Même si, arrivant également sur la dernière ligne droite, je n'ai pas beaucup de temps libre, j'y jetterai un oeil avec plaisir.