PGLS, model averaging fonction et sélection de modèle

Postez ici vos questions, réponses, commentaires ou suggestions - Les sujets seront ultérieurement répartis dans les archives par les modérateurs

Modérateur : Groupe des modérateurs

Théo Dokhelar
Messages : 2
Enregistré le : 10 Mar 2020, 08:05

PGLS, model averaging fonction et sélection de modèle

Messagepar Théo Dokhelar » 15 Juil 2020, 11:57

Bonjour,

Je travaille actuellement sur un jeu de données d'espèces d'oiseaux et cherche à voir sous R quels facteurs peuvent rendre certaines de ces espèces invasives.
Pour cela je travaille avec des PGLS qui sont des régressions linéaires prenant en compte la phylogénie des espèces utilisées dans cette régression.

Parmi mes 478 espèces étudiées, seulement 150 d'entres elles sont considérées comme invasives contre 328 non invasives. Pour éviter un biais dans mes analyses liés à cette différence de taille dans ces échantillons, j'ai donc réalisé un bootstrapping afin d'extraire aléatoirement 150 espèces parmi les non invasives et les ajouter dans ma PGLS aux 150 autres considérées comme invasives et ainsi faire tourner ma régression sur ce total de 300 espèces.

Après avoir fait tourner ce bootstrapping 1000 fois et stocké mes résultats de PGLS dans une liste, j'utilise la fonction de model averaging "model.avg" du package MuMIn sur la liste de mes résultats de PGLS afin d'obtenir un modèle final sur lequel je puisse interpréter mes résultats.

Cependant, à chaque fois que réalise ces étapes, je me retrouve avec un modèle final complètement différent, avec des p-values, coefficients et facteurs significatifs complètement différents. Je sais qu'il est normal que les modèles finaux ne soient pas tout à fait similaire du à la sélection aléatoire dans mon bootstrapping mais je ne comprends pas pourquoi les résultats varient autant d'un modèle à l'autre. Je ne sais donc pas sur quel modèle me baser pour interpréter mes résultats.

Auriez-vous une piste pour m'aider ?
Les variables que j'utilise dans mon modèle sont soit numériques soit factorielles.
J'ai essayé d'utiliser les corrélations Browniennes et de Pagel avec un lambda de 1 mais obtient le même problème avec ces deux méthodes.
J'ai bien entendu réalisé chaque étape séparément pour voir d'où pouvait venir mon problème mais je n'arrive toujours pas à le régler.

Voici le code de ma boucle :

Code : Tout sélectionner

intro<-subset(EU, Introduced == "1")
non_intro<-subset(EU, Introduced =="0")
phy<-read.newick(file="list_taxon_fin.txt")
ListPagel = list()
ListBro=list()
for (i in 1:1000){
  #create a sample of 150 non introduced birds and create a dataframe with the introduced species
  non_intro_sub=sample_n(non_intro,size = 150)
  EU3=bind_rows(intro,non_intro_sub)
  EU_final=data.frame(EU3,row.names=1)
 
  #creation fo the variables for the pgls
  introduced<-EU_final[,"Introduced"]
  diet<-EU_final[,"DietBreadth"]
  omni<-EU_final[,"Omnivorous"]
  carni<-EU_final[,"Carnivorous"]
  foli<-EU_final[,"Folivorous"]
  scav<-EU_final[,"Scavenger"]
  pisci<-EU_final[,"Piscivorous"]
  grani<-EU_final[,"Granivorous"]
  insecti<-EU_final[,"Insectivorous"]
  noc<-EU_final[,"Nocturnal"]
  color<-EU_final[,"Plumage"]
  status<-EU_final[,"Status"]
  niche<-EU_final[,"Niche_breadth"]
  hab<-EU_final[,"Habitat_position"]
  urban<-EU_final[,"Urban"]
  pet<-EU_final[,"Pet"]
  trade_nat<-EU_final[,"Trade..national."]
  hunt<-EU_final[,"Sport.hunting"]
  food<-EU_final[,"Food"]
  cites<-EU_final[,"Cites"]
  migr<-EU_final[,"migration"]
  area<-EU_final[,"Shape_Area"]
  temp<-EU_final[,"temprangemean"]
  HFP<-EU_final[,"HFPmean"]
  PC1<-EU_final[,"PC1"]
  PC2<-EU_final[,"PC2"]
 
  #removing tip not used in the tree
  birdtree<-name.check(phy,EU_final)
  tree_bootstrap<-drop.tip(phy,birdtree$tree_not_data)
 
  #correlation coef
  pagel_coef<-corPagel(1, tree_bootstrap)
  bro_coef<-corBrownian(1, tree_bootstrap)
 
  #model
  bootmodelPagel <- gls(introduced ~ diet+omni+carni+foli+scav+pisci+insecti+noc+color+status+niche+hab+urban+pet+trade_nat+hunt+food+cites+migr+area+temp+HFP+PC1+PC2, correlation = pagel_coef,
                   data = EU_final, method = "ML")
  bootmodelBro<- gls(introduced ~ diet+omni+carni+foli+scav+pisci+insecti+noc+color+status+niche+hab+urban+pet+trade_nat+hunt+food+cites+migr+area+temp+HFP+PC1+PC2, correlation = bro_coef,
                    data = EU_final, method = "ML")
  #to save the different models
  ListPagel[[length(ListPagel)+1]] = bootmodelPagel
  ListBro[[length(ListBro)+1]] = bootmodelBro
}
#model averaging
modavPagel=model.avg(ListPagel)
modavBro=model.avg(ListBro)

Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Google [Bot] et 1 invité