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)