####################################################################### ####################################################################### ##### ##### ##### ##### ##### Heritability Local Fitness Data Clownfish Kimbe Island ##### ##### <°)))>< ##### ##### ##### ####################################################################### ####################################################################### #################################### ############ PACKAGES ############## library(tnet) library(pedantics) library(ggplot2) library(wesanderson) library(scales) library(RColorBrewer) library(gplots) library(shape) library(proxy) library(MCMCglmm) library(nadiv) library(igraph) library(lme4) library(lmerTest) library(arm) library(sp) library(spdep) library(splm) library(ape) library(vegan) library(ade4) library(Cairo) library(geoR) library(proxy) library(pedigree) library(FactoMineR) library(MasterBayes) library(QGglmm) library(HDInterval) #################################################### ############ PEDIGREE + INVERSE #################### Ped=as.data.frame(read.table(file="D:/Documents/POSTDOC_TOULOUSE/Publications/TransmissionNiche/Data/PED_CLOWNFISH_FINAL_sansNR.txt",header=TRUE, na.string="NA")) Ped$ID=as.factor(Ped$ID) Ped$FATHER=as.factor(Ped$FATHER) Ped$MOTHER=as.factor(Ped$MOTHER) #Orders the pedigree so that parents come before their offspring library(MasterBayes) Ped=orderPed(Ped) #possibility to order the pedigree using fixPedigree if orderPed doesnt work... Ped=fixPedigree(Ped) for(x in 1:3) Ped[,x]<-as.factor(Ped[,x]) write.table(Ped, "D:/Documents/POSTDOC_TOULOUSE/Publications/TransmissionNiche/Data/OrderPed.txt", sep="\t") #we keep only informative individuals from the pedigree Ped=prunePed(Ped, Ped$ID) #Gets the inverse of the relatedness matrix library(pedantics) Ainv=inverseA(Ped)$Ainv drawPedigree(Ped) #draw the matrix of genetic similarity library(nadiv) pedigree=prepPed(Ped) #save the figure jpeg(filename="Models/Figs/Pedigreeall.jpg", width=900, height=1000,units="px") drawPedigree(Ped) dev.off() ################################################### ############ DATA REPEATED MEASURES ############### setwd("D:/Documents/POSTDOC_TOULOUSE/Publications/TransmissionNiche/Data") Aperc=as.data.frame(read.table(file="D:/Documents/POSTDOC_TOULOUSE/Publications/TransmissionNiche/Data/DATA_CLOWNFISH_FINAL_sansNR_REPEATED_MEASURES.txt",header=TRUE,na.string="NA")) #gets ID/dam/sire as factors Aperc$ANIMAL=as.factor(Aperc$ID) Aperc$MOTHER=as.factor(Aperc$MOTHER) Aperc$FATHER=as.factor(Aperc$FATHER) #these need to be repeated for connection to the extra matrices in the models Aperc$ANIMAL2=as.factor(Aperc$ID) Aperc$animal=as.factor(Aperc$ID) Aperc$MOTHER2=as.factor(Aperc$MOTHER) Aperc$FATHER2=as.factor(Aperc$FATHER) #these need to be repeated for connection to the extra matrices in the models Aperc$ANIMAL2=as.factor(Aperc$ID) Aperc$MOTHER2=as.factor(Aperc$MOTHER) Aperc$FATHER2=as.factor(Aperc$FATHER) Aperc$HABITATVIE=as.factor(Aperc$HABITATVIE) Aperc$HABITATMOTHER=as.factor(Aperc$HABITATMOTHER) # gets factor or numeric for each parameters Aperc$sex=as.factor(Aperc$sex) Aperc$cohort=as.factor(Aperc$cohort) Aperc$SuccessReproAR=as.numeric(Aperc$SuccessReproAR) Aperc$Survival=as.numeric(Aperc$Survival) Aperc$LRS=as.numeric(Aperc$LRS) Aperc$occurrence=as.numeric(Aperc$occurrence) Aperc$HabitatBirth=as.numeric(Aperc$HabitatBirth) Aperc$HabitatResident=as.numeric(Aperc$HabitatResident) Aperc$LRSbin=as.numeric(Aperc$LRS) Aperc$LRSbin[Aperc$LRSbin>0]<-1 Aperc$delifing=as.numeric(Aperc$delifing) Aperc$DepthAnemone=as.numeric(Aperc$DepthAnemone) Aperc$CodespAnemone=as.factor(Aperc$CodespAnemone) Aperc$CodeArea=as.factor(Aperc$CodeArea) Aperc$PropDeadCorals=as.numeric(Aperc$PropDeadCorals) Aperc$PropCorals=as.numeric(Aperc$PropCorals) Aperc$PropSand=as.numeric(Aperc$PropSand) Aperc$PropRock=as.numeric(Aperc$PropRock) Aperc$PropAlgae=as.numeric(Aperc$PropAlgae) Aperc$DepthAnemoneParents=as.numeric(Aperc$DepthAnemoneParents) Aperc$CodeAreaParents=as.factor(Aperc$CodeAreaParents) Aperc$CodespAnemoneParents=as.factor(Aperc$CodespAnemoneParents) Aperc$PropDeadCoralsParents=as.numeric(Aperc$PropDeadCoralsParents) Aperc$PropCoralsParents=as.numeric(Aperc$PropCoralsParents) Aperc$PropSandParents=as.numeric(Aperc$PropSandParents) Aperc$PropRockParents=as.numeric(Aperc$PropRockParents) Aperc$PropAlgaeParents=as.numeric(Aperc$PropAlgaeParents) #Visualisation des données + sauvegarde setwd("D:/Documents/POSTDOC_TOULOUSE/Publications/TransmissionNiche/Data") jpeg("Models/Figs/delifingRM.jpg") hist(Aperc$delifing) dev.off() jpeg("Models/Figs/histLRSRM.jpg") hist(Aperc$LRS) dev.off() jpeg("Models/Figs/LRSbinRM.jpg") hist(Aperc$LRSbin) dev.off() ################################### ############ MODELS ############### ##########################" #### PCI PAPER library(MCMCglmm) prior5=list(G=list(G1=list(V=1, nu=0.02), G2=list(V=1, nu=0.02), G3=list(V=1, nu=0.02), G4=list(V=1, nu=0.02), G5=list(V=1, nu=0.02)), R=list(V=1, nu=0.02)) MCMCglmm_LRS_pnas5<-MCMCglmm(LRS ~ 1 , random=~ANIMAL + MOTHER + HabitatBirth + HabitatResident + ID , family="poisson" , data= Aperc , ginverse=list(ANIMAL=Ainv) , prior=prior5 , nitt = 1000000 , burnin = 1000 , thin = 1000) save(MCMCglmm_LRS_pnas5, file=paste("Models/","Figs/","MCMCglmm_LRS_pnas5.Rdata",sep="")) prior4=list(G=list(G1=list(V=1, nu=0.02), G2=list(V=1, nu=0.02), G3=list(V=1, nu=0.02), G4=list(V=1, nu=0.02)), R=list(V=1, nu=0.02)) MCMCglmm_LRS_pnas4<-MCMCglmm(LRS ~ 1 , random=~ANIMAL + MOTHER + HabitatBirth + ID , family="poisson" , data= Aperc , ginverse=list(ANIMAL=Ainv) , prior=prior4 , nitt = 1000000 , burnin = 1000 , thin = 1000) save(MCMCglmm_LRS_pnas4, file=paste("Models/","Figs/","MCMCglmm_LRS_pnas4.Rdata",sep="")) MCMCglmm_LRS_pnas4b<-MCMCglmm(LRS ~ 1 , random=~ANIMAL + MOTHER + HabitatResident + ID , family="poisson" , data= Aperc , ginverse=list(ANIMAL=Ainv) , prior=prior4 , nitt = 1000000 , burnin = 1000 , thin = 1000) save(MCMCglmm_LRS_pnas4b, file=paste("Models/","Figs/","MCMCglmm_LRS_pnas4b.Rdata",sep="")) MCMCglmm_LRS_pnas4c<-MCMCglmm(LRS ~ 1 , random=~ANIMAL + HabitatBirth + HabitatResident + ID , family="poisson" , data= Aperc , ginverse=list(ANIMAL=Ainv) , prior=prior4 , nitt = 1000000 , burnin = 1000 , thin = 1000) save(MCMCglmm_LRS_pnas4c, file=paste("Models/","Figs/","MCMCglmm_LRS_pnas4c.Rdata",sep="")) MCMCglmm_LRS_pnas4d<-MCMCglmm(LRS ~ 1 , random=~MOTHER + HabitatBirth + HabitatResident + ID , family="poisson" , data= Aperc , ginverse=list(ANIMAL=Ainv) , prior=prior4 , nitt = 1000000 , burnin = 1000 , thin = 1000) save(MCMCglmm_LRS_pnas4d, file=paste("Models/","Figs/","MCMCglmm_LRS_pnas4d.Rdata",sep="")) ############## MCMCglmm_delifing_pnas5<-MCMCglmm(delifing ~ 1 , random=~ANIMAL + MOTHER + HabitatBirth + HabitatResident + ID , family="gaussian" , data= Aperc , ginverse=list(ANIMAL=Ainv) , prior=prior5 , nitt = 1000000 , burnin = 1000 , thin = 1000) save(MCMCglmm_delifing_pnas5, file=paste("Models/","Figs/","MCMCglmm_delifing_pnas5.Rdata",sep="")) prior4=list(G=list(G1=list(V=1, nu=0.02), G2=list(V=1, nu=0.02), G3=list(V=1, nu=0.02), G4=list(V=1, nu=0.02)), R=list(V=1, nu=0.02)) MCMCglmm_delifing_pnas4<-MCMCglmm(delifing ~ 1 , random=~ANIMAL + MOTHER + HabitatBirth + ID , family="gaussian" , data= Aperc , ginverse=list(ANIMAL=Ainv) , prior=prior4 , nitt = 1000000 , burnin = 1000 , thin = 1000) save(MCMCglmm_delifing_pnas4, file=paste("Models/","Figs/","MCMCglmm_delifing_pnas4.Rdata",sep="")) MCMCglmm_delifing_pnas4b<-MCMCglmm(delifing ~ 1 , random=~ANIMAL + MOTHER + HabitatResident + ID , family="gaussian" , data= Aperc , ginverse=list(ANIMAL=Ainv) , prior=prior4 , nitt = 1000000 , burnin = 1000 , thin = 1000) save(MCMCglmm_delifing_pnas4b, file=paste("Models/","Figs/","MCMCglmm_delifing_pnas4b.Rdata",sep="")) MCMCglmm_delifing_pnas4c<-MCMCglmm(delifing ~ 1 , random=~ANIMAL + HabitatBirth + HabitatResident + ID , family="gaussian" , data= Aperc , ginverse=list(ANIMAL=Ainv) , prior=prior4 , nitt = 1000000 , burnin = 1000 , thin = 1000) save(MCMCglmm_delifing_pnas4c, file=paste("Models/","Figs/","MCMCglmm_delifing_pnas4c.Rdata",sep="")) MCMCglmm_delifing_pnas4d<-MCMCglmm(delifing ~ 1 , random=~MOTHER + HabitatBirth + HabitatResident + ID , family="gaussian" , data= Aperc , ginverse=list(ANIMAL=Ainv) , prior=prior4 , nitt = 1000000 , burnin = 1000 , thin = 1000) save(MCMCglmm_delifing_pnas4d, file=paste("Models/","Figs/","MCMCglmm_delifing_pnas4d.Rdata",sep="")) setwd("D:/Documents/POSTDOC_TOULOUSE/Publications/TransmissionNiche/Data/Models/Figs") load("MCMCglmm_delifing_pnas5.Rdata") load("MCMCglmm_delifing_pnas4.Rdata") load("MCMCglmm_delifing_pnas4.Rdata") load("MCMCglmm_delifing_pnas4b.Rdata") load("MCMCglmm_delifing_pnas4c.Rdata") load("MCMCglmm_delifing_pnas4d.Rdata") load("MCMCglmm_LRS_pnas5.Rdata") load("MCMCglmm_LRS_pnas4.Rdata") load("MCMCglmm_LRS_pnas4.Rdata") load("MCMCglmm_LRS_pnas4b.Rdata") load("MCMCglmm_LRS_pnas4c.Rdata") load("MCMCglmm_LRS_pnas4d.Rdata") summary(MCMCglmm_delifing_pnas5) posterior.mode(MCMCglmm_delifing_pnas5$VCV) HPDinterval(MCMCglmm_delifing_pnas5$VCV) summary(MCMCglmm_delifing_pnas4) posterior.mode(MCMCglmm_delifing_pnas4$VCV) HPDinterval(MCMCglmm_delifing_pnas4$VCV) VA_delifing4= MCMCglmm_delifing_pnas4$VCV[,"ANIMAL"] VM_delifing4=MCMCglmm_delifing_pnas4$VCV[,"MOTHER"] VNH_delifing4=MCMCglmm_delifing_pnas4$VCV[,"HabitatBirth"] VPE_delifing4=MCMCglmm_delifing_pnas4$VCV[,"ID"] VR_delifing4=MCMCglmm_delifing_pnas4$VCV[,"units"] #heritability H2_delifing4=VA_delifing4/(VA_delifing4 + VM_delifing4 + VNH_delifing4 + VPE_delifing4 + VR_delifing4) posterior.mode(H2_delifing4) HPDinterval(H2_delifing4,0.95) summary(MCMCglmm_delifing_pnas4) summary(MCMCglmm_delifing_pnas4b) posterior.mode(MCMCglmm_delifing_pnas4b$VCV) HPDinterval(MCMCglmm_delifing_pnas4b$VCV) VA_delifing4b= MCMCglmm_delifing_pnas4b$VCV[,"ANIMAL"] VM_delifing4b=MCMCglmm_delifing_pnas4b$VCV[,"MOTHER"] VNH_delifing4b=MCMCglmm_delifing_pnas4b$VCV[,"HabitatBirth"] VRH_delifing4b=MCMCglmm_delifing_pnas4b$VCV[,"HabitatResident"] VPE_delifing4b=MCMCglmm_delifing_pnas4b$VCV[,"ID"] VR_delifing4b=MCMCglmm_delifing_pnas4b$VCV[,"units"] #heritability H2_delifing4b=VA_delifing4b/(VA_delifing4b + VM_delifing4b + VRH_delifing4b + VPE_delifing4b + VR_delifing4b) posterior.mode(H2_delifing4b) HPDinterval(H2_delifing4b,0.95) summary(MCMCglmm_delifing_pnas4c) posterior.mode(MCMCglmm_delifing_pnas4c$VCV) HPDinterval(MCMCglmm_delifing_pnas4c$VCV) VA_delifing4c= MCMCglmm_delifing_pnas4c$VCV[,"ANIMAL"] VM_delifing4c=MCMCglmm_delifing_pnas4c$VCV[,"MOTHER"] VNH_delifing4c=MCMCglmm_delifing_pnas4c$VCV[,"HabitatBirth"] VRH_delifing4c=MCMCglmm_delifing_pnas4c$VCV[,"HabitatResident"] VPE_delifing4c=MCMCglmm_delifing_pnas4c$VCV[,"ID"] VR_delifing4c=MCMCglmm_delifing_pnas4c$VCV[,"units"] #heritability H2_delifing4c=VA_delifing4c/(VA_delifing4c + VNH_delifing4c + VRH_delifing4c+ VPE_delifing4c + VR_delifing4c) posterior.mode(H2_delifing4c) HPDinterval(H2_delifing4c,0.95) summary(MCMCglmm_delifing_pnas4d) posterior.mode(MCMCglmm_delifing_pnas4d$VCV) HPDinterval(MCMCglmm_delifing_pnas4d$VCV) summary(MCMCglmm_LRS_pnas5) posterior.mode(MCMCglmm_LRS_pnas5$VCV) HPDinterval(MCMCglmm_LRS_pnas5$VCV) summary(MCMCglmm_LRS_pnas4) posterior.mode(MCMCglmm_LRS_pnas4$VCV) HPDinterval(MCMCglmm_LRS_pnas4$VCV) summary(MCMCglmm_LRS_pnas4b) posterior.mode(MCMCglmm_LRS_pnas4b$VCV) HPDinterval(MCMCglmm_LRS_pnas4b$VCV) summary(MCMCglmm_LRS_pnas4c) posterior.mode(MCMCglmm_LRS_pnas4c$VCV) HPDinterval(MCMCglmm_LRS_pnas4c$VCV) summary(MCMCglmm_LRS_pnas4d) posterior.mode(MCMCglmm_LRS_pnas4d$VCV) HPDinterval(MCMCglmm_LRS_pnas4d$VCV) load("MCMCglmm_delifing_pnas5.Rdata") summary(MCMCglmm_delifing_pnas5) VA_delifing= MCMCglmm_delifing_pnas5$VCV[,"ANIMAL"] mean(VA_delifing) VM_delifing=MCMCglmm_delifing_pnas5$VCV[,"MOTHER"] mean(VM_delifing) VNH_delifing=MCMCglmm_delifing_pnas5$VCV[,"HabitatBirth"] mean(VNH_delifing) VRH_delifing=MCMCglmm_delifing_pnas5$VCV[,"HabitatResident"] mean(VRH_delifing) VPE_delifing=MCMCglmm_delifing_pnas5$VCV[,"ID"] mean(VPE_delifing) VR_delifing=MCMCglmm_delifing_pnas5$VCV[,"units"] mean(VR_delifing) posterior.mode(MCMCglmm_delifing_pnas5$VCV) HPDinterval(MCMCglmm_delifing_pnas5$VCV) mean(MCMCglmm_delifing_pnas5$VCV) #heritability H2_delifing=VA_delifing/(VA_delifing + VM_delifing + VNH_delifing + VRH_delifing + VPE_delifing + VR_delifing) mean(H2_delifing) posterior.mode(H2_delifing) HPDinterval(H2_delifing,0.95) #mother effects M2_delifing=VM_delifing/(VA_delifing + VM_delifing + VNH_delifing + VRH_delifing + VPE_delifing + VR_delifing) posterior.mode(M2_delifing) HPDinterval(M2_delifing,0.95) mean(M2_delifing) #evolvability mdelifing=mean(Aperc$delifing) Ia_delifing=VA_delifing/(mdelifing*mdelifing) posterior.mode(Ia_delifing) HPDinterval(Ia_delifing,0.95) QGparams(mu=mu,var.a=va,var.p=vp,model="Poisson.log") #regarder si mean.obs et var.obs sont proches des données obs mean(Aperc$LRS) var(Aperc$LRS) summary(MCMCglmm_delifing_pnas5) # model complets MCMCglmm_LRS_pnas5 et MCMCglmm_delifing_pnas5 #Diagnostic of the MCMC setwd("D:/Documents/POSTDOC_TOULOUSE/Publications/TransmissionNiche/Data") jpeg("Models/Figs/MCMCglmm_LRS_PLOTVCV.jpg") plot(MCMCglmm_LRS_pnas5$VCV) dev.off() jpeg("Models/Figs/MCMCglmm_delifing_PLOTVCV.jpg") plot(MCMCglmm_delifing_pnas5$VCV) dev.off() #Results from the MCMCglmm posterior.mode(MCMCglmm_LRS_pnas5$VCV) HPDinterval(MCMCglmm_LRS_pnas5$VCV) summary(MCMCglmm_LRS_pnas5) VA_LRS= MCMCglmm_LRS_pnas5$VCV[,"ANIMAL"] VM_LRS=MCMCglmm_LRS_pnas5$VCV[,"MOTHER"] VNH_LRS=MCMCglmm_LRS_pnas5$VCV[,"HabitatBirth"] VRH_LRS=MCMCglmm_LRS_pnas5$VCV[,"HabitatResident"] VPE_LRS=MCMCglmm_LRS_pnas5$VCV[,"ID"] VR_LRS=MCMCglmm_LRS_pnas5$VCV[,"units"] mean(VA_LRS) mean(VM_LRS) mean(VNH_LRS) mean(VRH_LRS) mean(VPE_LRS) mean(VR_LRS) #heritability H2_LRS=VA_LRS/(VA_LRS + VM_LRS + VNH_LRS + VRH_LRS + VPE_LRS + VR_LRS) mean(H2_LRS) posterior.mode(H2_LRS) HPDinterval(H2_LRS,0.95) #mother effects M2_LRS=VM_LRS/(VA_LRS + VM_LRS + VNH_LRS + VRH_LRS + VPE_LRS + VR_LRS) mean(M2_LRS) posterior.mode(M2_LRS) HPDinterval(M2_LRS,0.95) #evolvability equation: Ia=Va/m² with m the mean of the trait mLRS=mean(Aperc$LRS) Ia_LRS=mean(VA_LRS)/(mLRS*mLRS) mean(Ia_LRS) posterior.mode(Ia_LRS) HPDinterval(Ia_LRS,0.95) library(QGglmm) library(MCMCglmm) #heritability estimate using de Villmereuille et al 2016 Genetics mu= mean(MCMCglmm_LRS_pnas5[["Sol"]][,"(Intercept)"]) va= mean(MCMCglmm_LRS_pnas5[["VCV"]][,"ANIMAL"]) vm=mean(MCMCglmm_LRS_pnas5[["VCV"]][,"MOTHER"]) vnh=mean(MCMCglmm_LRS_pnas5[["VCV"]][,"HabitatBirth"]) vrh=mean(MCMCglmm_LRS_pnas5[["VCV"]][,"HabitatResident"]) vpe=mean(MCMCglmm_LRS_pnas5[["VCV"]][,"ID"]) vr=mean(MCMCglmm_LRS_pnas5[["VCV"]][,"units"]) vp=mean(rowSums(MCMCglmm_LRS_pnas5[["VCV"]])) QGparams(mu=mu,var.a=va,var.p=vp,model="Poisson.log") #héritabilité QGparams(mu=mu,var.a=vm,var.p=vp,model="Poisson.log") # effets maternels QGparams(mu=mu,var.a=vnh,var.p=vp,model="Poisson.log") #HabitatBirth QGparams(mu=mu,var.a=vrh,var.p=vp,model="Poisson.log") #HabitatResident QGparams(mu=mu,var.a=vpe,var.p=vp,model="Poisson.log") #ID QGparams(mu=mu,var.a=vr,var.p=vp,model="Poisson.log") #units QGicc(mu=mu,var.comp=vm,var.p=vp,model="Poisson.log") # effets maternels QGicc(mu=mu,var.comp=vnh,var.p=vp,model="Poisson.log") #HabitatBirth QGicc(mu=mu,var.comp=vrh,var.p=vp,model="Poisson.log") #HabitatResident QGicc(mu=mu,var.comp=vpe,var.p=vp,model="Poisson.log") #ID QGicc(mu=mu,var.comp=vr,var.p=vp,model="Poisson.log") #units QGicc(mu = mu, var.p = vp, var.comp = va, model = "Poisson.log")#Computing the broad - sense heritability QGicc(mu = mu, var.p = vp, var.comp = vm, model = "Poisson.log") #réccupérer les posteriors distributions df= data.frame(mu = as.vector(MCMCglmm_LRS_pnas5[["Sol"]][, "(Intercept)"]), va = as.vector(MCMCglmm_LRS_pnas5[["VCV"]][, "ANIMAL"]), vm= as.vector(MCMCglmm_LRS_pnas5[["VCV"]][, "MOTHER"]), vnh=mean(MCMCglmm_LRS_pnas5[["VCV"]][,"HabitatBirth"]), vrh=mean(MCMCglmm_LRS_pnas5[["VCV"]][,"HabitatResident"]), vpe=mean(MCMCglmm_LRS_pnas5[["VCV"]][,"ID"]), vr=mean(MCMCglmm_LRS_pnas5[["VCV"]][,"units"]), vp = rowSums(MCMCglmm_LRS_pnas5[["VCV"]])) head(df) # on calcule les posterior distribution pour h² post= do.call("rbind", apply(df, 1, function(row){ QGparams(mu=row[["mu"]], var.a=row[["va"]], var.p=row[["vp"]], model="Poisson.log", verbose= FALSE) })) head(post) # on regarde la tête du tableau post.h2.obs=as.matrix(post$h2.obs) # on sélectionne la colonne h² mean(post.h2.obs) post.var.a.obs=as.matrix(post$var.a.obs) hdi(post.var.a.obs, credMass = 0.95)# on mesure la CI95% hdi(post.h2.obs, credMass = 0.95) #estimation de l'évolvabilité du trait LRS post_Ia_LRS_obervedscale= post$var.a.obs / (mLRS*mLRS) mean(post_Ia_LRS_obervedscale) posterior.mode(post_Ia_LRS_obervedscale) hdi(post_Ia_LRS_obervedscale,0.95) # on calcule les posterior distribution pour m² postm= do.call("rbind", apply(df, 1, function(row){ QGparams(mu=row[["mu"]], var.a=row[["vm"]], var.p=row[["vp"]], model="Poisson.log", verbose= FALSE) })) head(postm) # on regarde la tête du tableau postm.h2.obs=as.matrix(postm$h2.obs) # on sélectionne la colonne h² postm.var.a.obs=as.matrix(postm$var.a.obs) hdi(postm.h2.obs, credMass = 0.95)# on mesure la CI95% hdi(postm.var.a.obs, credMass = 0.95) # on calcule les posterior distribution pour m² avec QGicc postm= do.call("rbind", apply(df, 1, function(row){ QGicc(mu=row[["mu"]], var.comp=row[["vm"]], var.p=row[["vp"]], model="Poisson.log", verbose= FALSE) })) head(postm) # on regarde la tête du tableau postm.var.comp.obs=as.matrix(postm$var.comp.obs) hdi(postm.var.comp.obs, credMass = 0.95) postm.m2.obs=as.matrix(postm$icc.obs) # on sélectionne la colonne h² hdi(postm.m2.obs, credMass = 0.95)# on mesure la CI95% # on calcule les posterior distribution pour vnh postvnh= do.call("rbind", apply(df, 1, function(row){ QGicc(mu=row[["mu"]], var.comp=row[["vnh"]], var.p=row[["vp"]], model="Poisson.log", verbose= FALSE) })) head(postvnh) # on regarde la tête du tableau postvnh.var.comp.obs=as.matrix(postvnh$var.comp.obs) # on sélectionne la colonne h² hdi(postvnh.var.comp.obs, credMass = 0.95)# on mesure la CI95% # on calcule les posterior distribution pour vrh postvrh= do.call("rbind", apply(df, 1, function(row){ QGicc(mu=row[["mu"]], var.comp=row[["vrh"]], var.p=row[["vp"]], model="Poisson.log", verbose= FALSE) })) head(postvrh) # on regarde la tête du tableau postvrh.var.comp.obs=as.matrix(postvrh$var.comp.obs) # on sélectionne la colonne h² hdi(postvrh.var.comp.obs, credMass = 0.95)# on mesure la CI95% # on calcule les posterior distribution pour vpe postvpe= do.call("rbind", apply(df, 1, function(row){ QGicc(mu=row[["mu"]], var.comp=row[["vpe"]], var.p=row[["vp"]], model="Poisson.log", verbose= FALSE) })) head(postvpe) # on regarde la tête du tableau postvpe.var.comp.obs=as.matrix(postvpe$var.comp.obs) # on sélectionne la colonne h² hdi(postvpe.var.comp.obs, credMass = 0.95)# on mesure la CI95% # on calcule les posterior distribution pour vr postvr= do.call("rbind", apply(df, 1, function(row){ QGicc(mu=row[["mu"]], var.comp=row[["vr"]], var.p=row[["vp"]], model="Poisson.log", verbose= FALSE) })) head(postvr) # on regarde la tête du tableau postvr.var.comp.obs=as.matrix(postvr$var.comp.obs) # on sélectionne la colonne h² hdi(postvr.var.comp.obs, credMass = 0.95)# on mesure la CI95% posterior.mode(MCMCglmm_LRS_pnas5$VCV)#estimation of additive genetics ans residual variance HPDinterval(MCMCglmm_LRS_pnas5$VCV) #interval de confiance (95%) du mod distribution. ######### Utiliser QGparams pour vérifier que Poisson log est bien pour LRS # QGparams renvoie moyenne et variance des données qui correspondent à l'estimation du modèle setwd("D:/Documents/POSTDOC_TOULOUSE/Publications/TransmissionNiche/Data/Models/Figs") load("MCMCglmm_LRS_pnas5.Rdata") QGparams(mu=mu,var.a=va,var.p=vp,model="Poisson.log") #regarder si mean.obs et var.obs sont proches des données obs mean(Aperc$LRS) var(Aperc$LRS) summary(MCMCglmm_delifing_pnas5) setwd("D:/Documents/POSTDOC_TOULOUSE/Publications/TransmissionNiche/Data/Models/Figs") ##################################################################################