####### lecture fichier donnees # ci-dessous : lecture via commandes R ... ne pas considerer les deux lignes suivantes si utilisation de "Import Dataset" de R-Studio #setwd("~/Documents/Enseignement/Master Agro/Cours") #Mais <- read.csv("Exemple_Mais.txt", sep="") options(contrasts = c("contr.treatment", "contr.poly")) # (par default) #options(contrasts = c("contr.sum", "contr.poly")) Mais <- read.csv("~/Documents/Enseignement/Master Agro/Cours/Exemple_Mais.txt", sep="") # jeu de donnees : Mais head(Mais) str(Mais) # transformation en facteurs Mais$Engrais<-as.factor(Mais$Engrais) Mais$Varietes<-as.factor(Mais$Varietes) xtabs(~Mais$Engrais+Mais$Variete) # graphes des reponses moyennes with(Mais, interaction.plot(Engrais,Varietes,prod)) # potentiellement de l'interaction ! #table anova ... plusieurs facons de l'obtenir selon l appel des vecteurs ! # summary(aov(Mais[,1]~Mais[,2]*Mais[,3])) #ou summary(aov(Mais$prod~Mais$Engrais*Mais$Varietes)) #ou summary(aov(prod~Engrais*Varietes,data=Mais)) #ou with(Mais, summary(aov(prod~Engrais*Varietes))) #ou mod_mais<-lm(prod~Engrais*Varietes,data=Mais) # ici modelisation ! #mod_mais<-lm(prod~Engrais+Varietes+Engrais:Varietes,data=Mais) anova(mod_mais) ######## # renomer les niveaux des facteurs pour plus de lisibilite des resultats # Engrais ---> E1 E2 E3 plutot que 1 2 3 # effectifs table(Mais$Engrais) # 1 2 3 # 8 8 8 # structure : chaque modalite est successivement repetee 8 fois dans le tableau de depart Mais$Engrais<-gl(3, 8,labels=c("E1","E2", "E3") ) # 3 labels repliques chacun 8 fois de suite Mais$Engrais # [1] E1 E1 E1 E1 E1 E1 E1 E1 E2 E2 E2 E2 E2 E2 E2 E2 E3 E3 E3 E3 E3 E3 E3 E3 # Levels: E1 E2 E3 # Varietes ---> V1 V2 V3 V4 plutot que 1 2 3 4 table(Mais$Varietes) #1 2 3 4 #6 6 6 6 # structure : ici chaque modalite est successivementrepetee 2 fois dans le tableau de depart Mais$Varietes<-gl(4,2,labels=c("V1","V2", "V3", "V4")) # 4 labels repliques chacun 2 fois de suite Mais$Varietes #[1] V1 V1 V2 V2 V3 V3 V4 V4 V1 V1 V2 V2 V3 V3 V4 V4 V1 V1 V2 V2 V3 V3 V4 V4 #Levels: V1 V2 V3 V4 # intercation significative donc on doit travailler # sur effets d'un facteur a niveau fixe de l'autre # visualiser reponses moyennes en fonction des traitements consideres # creer un vecteur "traitement" (ie croisement des modalites des deux facteurs) interac<-interaction(Mais$Engrais, Mais$Varietes, sep="_") interac #[1] E1_V1 E1_V1 E1_V1 E1_V1 E1_V1 E1_V1 E1_V2 E1_V2 E2_V2 E2_V2 E2_V2 E2_V2 E2_V3 E2_V3 E2_V3 # [16] E2_V3 E3_V3 E3_V3 E3_V4 E3_V4 E3_V4 E3_V4 E3_V4 E3_V4 # Levels: E1_V1 E2_V1 E3_V1 E1_V2 E2_V2 E3_V2 E1_V3 E2_V3 E3_V3 E1_V4 E2_V4 E3_V4 # on a bien 3*4=12 traitements (croisements des modalites Varietes et Engrais) # boxplots reponse selon modalités du facteur interaction with(Mais, boxplot(prod~interac)) # ---> on constate des differences potentielles dans les reponses moyennes selon les traitements # recherche effets principaux de l'un des facteurs conditionnellement à modialte fixe de l'autre #par ex : effet de Variete conditionnellement a modalite 1 de l'Engrais ? # Approche 1 : se ramener à une ANOVA un facteur avec facteur construit via les traitements #et utiliser les contrastes pour tester les relations d'intérêt # ---> creer un "nouveau" facteur traitement correspondant aux # croisements des modalites des deux facteurs et faire une anova # un seul facteur eu utilisant ce facteur # nouveau facteur est le facteur "interac" créé précédemment modINTER<-lm(Mais$prod~interac) anova(modINTER) # ---> ici anova a un seul facteur ! # on retrouve la meme p-value que celle du model complet avec interaction # pour verifier : summary(lm(Mais$prod~Mais$Engrais*Mais$Varietes)) # recherche de differences potentielles via des approches comparaisons multiples # des modalités du facteur interac pour reperer les effets d'un facteur conditionnelement # à niveaux de l'autre library(multcomp) mc_tukey <- glht(modINTER, linfct=mcp(interac="Tukey")) summary(mc_tukey) #Linear Hypotheses: # Estimate Std. Error t value Pr(>|t|) #E2_V1 - E1_V1 == 0 3.000e+00 9.789e-01 3.065 0.1980 #E3_V1 - E1_V1 == 0 3.000e+00 9.789e-01 3.065 0.1982 #E1_V2 - E1_V1 == 0 3.000e+00 9.789e-01 3.065 0.1991 #E2_V2 - E1_V1 == 0 9.000e+00 9.789e-01 9.194 <0.01 *** #E3_V2 - E1_V1 == 0 6.000e+00 9.789e-01 6.129 <0.01 ** #E1_V3 - E1_V1 == 0 -2.000e+00 9.789e-01 -2.043 0.6640 #E2_V3 - E1_V1 == 0 1.032e-15 9.789e-01 0.000 1.0000 #E3_V3 - E1_V1 == 0 3.500e+00 9.789e-01 3.575 0.0930 . #E1_V4 - E1_V1 == 0 0.000e+00 9.789e-01 0.000 1.0000 #E2_V4 - E1_V1 == 0 -3.000e+00 9.789e-01 -3.065 0.1983 #E3_V4 - E1_V1 == 0 -2.000e+00 9.789e-01 -2.043 0.6636 #E3_V1 - E2_V1 == 0 3.997e-15 9.789e-01 0.000 1.0000 #E1_V2 - E2_V1 == 0 -4.441e-16 9.789e-01 0.000 1.0000 #E2_V2 - E2_V1 == 0 6.000e+00 9.789e-01 6.129 <0.01 ** #E3_V2 - E2_V1 == 0 3.000e+00 9.789e-01 3.065 0.1987 #E1_V3 - E2_V1 == 0 -5.000e+00 9.789e-01 -5.108 <0.01 ** #E2_V3 - E2_V1 == 0 -3.000e+00 9.789e-01 -3.065 0.1970 #E3_V3 - E2_V1 == 0 5.000e-01 9.789e-01 0.511 1.0000 #E1_V4 - E2_V1 == 0 -3.000e+00 9.789e-01 -3.065 0.1988 #E2_V4 - E2_V1 == 0 -6.000e+00 9.789e-01 -6.129 <0.01 ** #E3_V4 - E2_V1 == 0 -5.000e+00 9.789e-01 -5.108 <0.01 ** #E1_V2 - E3_V1 == 0 -4.441e-15 9.789e-01 0.000 1.0000 #E2_V2 - E3_V1 == 0 6.000e+00 9.789e-01 6.129 <0.01 ** #E3_V2 - E3_V1 == 0 3.000e+00 9.789e-01 3.065 0.1979 #E1_V3 - E3_V1 == 0 -5.000e+00 9.789e-01 -5.108 <0.01 ** #E2_V3 - E3_V1 == 0 -3.000e+00 9.789e-01 -3.065 0.1989 #E3_V3 - E3_V1 == 0 5.000e-01 9.789e-01 0.511 1.0000 #E1_V4 - E3_V1 == 0 -3.000e+00 9.789e-01 -3.065 0.1990 #E2_V4 - E3_V1 == 0 -6.000e+00 9.789e-01 -6.129 <0.01 ** #E3_V4 - E3_V1 == 0 -5.000e+00 9.789e-01 -5.108 <0.01 ** #E2_V2 - E1_V2 == 0 6.000e+00 9.789e-01 6.129 <0.01 ** #E3_V2 - E1_V2 == 0 3.000e+00 9.789e-01 3.065 0.1986 #E1_V3 - E1_V2 == 0 -5.000e+00 9.789e-01 -5.108 <0.01 ** #E2_V3 - E1_V2 == 0 -3.000e+00 9.789e-01 -3.065 0.1988 #E3_V3 - E1_V2 == 0 5.000e-01 9.789e-01 0.511 1.0000 #E1_V4 - E1_V2 == 0 -3.000e+00 9.789e-01 -3.065 0.1992 #E2_V4 - E1_V2 == 0 -6.000e+00 9.789e-01 -6.129 <0.01 ** #E3_V4 - E1_V2 == 0 -5.000e+00 9.789e-01 -5.108 <0.01 ** #E3_V2 - E2_V2 == 0 -3.000e+00 9.789e-01 -3.065 0.1981 #E1_V3 - E2_V2 == 0 -1.100e+01 9.789e-01 -11.237 <0.01 *** #E2_V3 - E2_V2 == 0 -9.000e+00 9.789e-01 -9.194 <0.01 *** #E3_V3 - E2_V2 == 0 -5.500e+00 9.789e-01 -5.618 <0.01 ** #E1_V4 - E2_V2 == 0 -9.000e+00 9.789e-01 -9.194 <0.01 *** #E2_V4 - E2_V2 == 0 -1.200e+01 9.789e-01 -12.258 <0.01 *** #E3_V4 - E2_V2 == 0 -1.100e+01 9.789e-01 -11.237 <0.01 *** #E1_V3 - E3_V2 == 0 -8.000e+00 9.789e-01 -8.172 <0.01 *** #E2_V3 - E3_V2 == 0 -6.000e+00 9.789e-01 -6.129 <0.01 ** #E3_V3 - E3_V2 == 0 -2.500e+00 9.789e-01 -2.554 0.3908 #E1_V4 - E3_V2 == 0 -6.000e+00 9.789e-01 -6.129 <0.01 ** #E2_V4 - E3_V2 == 0 -9.000e+00 9.789e-01 -9.194 <0.01 *** #E3_V4 - E3_V2 == 0 -8.000e+00 9.789e-01 -8.172 <0.01 *** #E2_V3 - E1_V3 == 0 2.000e+00 9.789e-01 2.043 0.6638 #E3_V3 - E1_V3 == 0 5.500e+00 9.789e-01 5.618 <0.01 ** #E1_V4 - E1_V3 == 0 2.000e+00 9.789e-01 2.043 0.6638 #E2_V4 - E1_V3 == 0 -1.000e+00 9.789e-01 -1.022 0.9935 #E3_V4 - E1_V3 == 0 -7.105e-15 9.789e-01 0.000 1.0000 #E3_V3 - E2_V3 == 0 3.500e+00 9.789e-01 3.575 0.0923 . #E1_V4 - E2_V3 == 0 -1.032e-15 9.789e-01 0.000 1.0000 #E2_V4 - E2_V3 == 0 -3.000e+00 9.789e-01 -3.065 0.1981 #E3_V4 - E2_V3 == 0 -2.000e+00 9.789e-01 -2.043 0.6641 #E1_V4 - E3_V3 == 0 -3.500e+00 9.789e-01 -3.575 0.0930 . #E2_V4 - E3_V3 == 0 -6.500e+00 9.789e-01 -6.640 <0.01 *** #E3_V4 - E3_V3 == 0 -5.500e+00 9.789e-01 -5.618 <0.01 ** #E2_V4 - E1_V4 == 0 -3.000e+00 9.789e-01 -3.065 0.1986 #E3_V4 - E1_V4 == 0 -2.000e+00 9.789e-01 -2.043 0.6640 #E3_V4 - E2_V4 == 0 1.000e+00 9.789e-01 1.022 0.9934 #--- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #(Adjusted p values reported -- single-step method) par(mar=c(3,7,3,3)) # par defaut : c(bottom, left, top, right)=c(5, 4, 4, 2) + 0.1 plot(mc_tukey) # approche SNK sur le facteur interac comparSNK<-agricolae::SNK.test(modINTER,"interac") comparSNK #$groups #Mais$prod groups #E2_V2 23.5 a #E3_V2 20.5 b #E3_V3 18.0 c #E1_V2 17.5 c #E2_V1 17.5 c #E3_V1 17.5 c #E1_V1 14.5 d #E1_V4 14.5 d #E2_V3 14.5 d #E1_V3 12.5 d #E3_V4 12.5 d #E2_V4 11.5 d # on pourrait aussi utiliser la fonction TukeyHSD deja vue mais mais ne s'applique que # sur objet aov (pas lm) donc obligation de refaire tourner aov() ... # en fait, idem approche Tukey avec library multcomp modINTER2<-aov(Mais$prod~interac) TukeyHSD(modINTER2) #$interac #diff lwr upr p adj #E2_V1-E1_V1 3.000000e+00 -0.8865633 6.8865633 0.1985174 #E3_V1-E1_V1 3.000000e+00 -0.8865633 6.8865633 0.1985174 #E1_V2-E1_V1 3.000000e+00 -0.8865633 6.8865633 0.1985174 #E2_V2-E1_V1 9.000000e+00 5.1134367 12.8865633 0.0000350 #E3_V2-E1_V1 6.000000e+00 2.1134367 9.8865633 0.0018056 #E1_V3-E1_V1 -2.000000e+00 -5.8865633 1.8865633 0.6637770 #E2_V3-E1_V1 -1.776357e-15 -3.8865633 3.8865633 1.0000000 #E3_V3-E1_V1 3.500000e+00 -0.3865633 7.3865633 0.0926317 #E1_V4-E1_V1 -3.552714e-15 -3.8865633 3.8865633 1.0000000 #E2_V4-E1_V1 -3.000000e+00 -6.8865633 0.8865633 0.1985174 #E3_V4-E1_V1 -2.000000e+00 -5.8865633 1.8865633 0.6637770 #E3_V1-E2_V1 3.552714e-15 -3.8865633 3.8865633 1.0000000 #E1_V2-E2_V1 3.552714e-15 -3.8865633 3.8865633 1.0000000 #E2_V2-E2_V1 6.000000e+00 2.1134367 9.8865633 0.0018056 #E3_V2-E2_V1 3.000000e+00 -0.8865633 6.8865633 0.1985174 #E1_V3-E2_V1 -5.000000e+00 -8.8865633 -1.1134367 0.0083586 #E2_V3-E2_V1 -3.000000e+00 -6.8865633 0.8865633 0.1985174 #E3_V3-E2_V1 5.000000e-01 -3.3865633 4.3865633 0.9999877 #E1_V4-E2_V1 -3.000000e+00 -6.8865633 0.8865633 0.1985174 #E2_V4-E2_V1 -6.000000e+00 -9.8865633 -2.1134367 0.0018056 #E3_V4-E2_V1 -5.000000e+00 -8.8865633 -1.1134367 0.0083586 #E1_V2-E3_V1 0.000000e+00 -3.8865633 3.8865633 1.0000000 #E2_V2-E3_V1 6.000000e+00 2.1134367 9.8865633 0.0018056 #E3_V2-E3_V1 3.000000e+00 -0.8865633 6.8865633 0.1985174 #E1_V3-E3_V1 -5.000000e+00 -8.8865633 -1.1134367 0.0083586 #E2_V3-E3_V1 -3.000000e+00 -6.8865633 0.8865633 0.1985174 #E3_V3-E3_V1 5.000000e-01 -3.3865633 4.3865633 0.9999877 #E1_V4-E3_V1 -3.000000e+00 -6.8865633 0.8865633 0.1985174 #E2_V4-E3_V1 -6.000000e+00 -9.8865633 -2.1134367 0.0018056 #E3_V4-E3_V1 -5.000000e+00 -8.8865633 -1.1134367 0.0083586 #E2_V2-E1_V2 6.000000e+00 2.1134367 9.8865633 0.0018056 #E3_V2-E1_V2 3.000000e+00 -0.8865633 6.8865633 0.1985174 #E1_V3-E1_V2 -5.000000e+00 -8.8865633 -1.1134367 0.0083586 #E2_V3-E1_V2 -3.000000e+00 -6.8865633 0.8865633 0.1985174 #E3_V3-E1_V2 5.000000e-01 -3.3865633 4.3865633 0.9999877 #E1_V4-E1_V2 -3.000000e+00 -6.8865633 0.8865633 0.1985174 #E2_V4-E1_V2 -6.000000e+00 -9.8865633 -2.1134367 0.0018056 #E3_V4-E1_V2 -5.000000e+00 -8.8865633 -1.1134367 0.0083586 #E3_V2-E2_V2 -3.000000e+00 -6.8865633 0.8865633 0.1985174 #E1_V3-E2_V2 -1.100000e+01 -14.8865633 -7.1134367 0.0000041 #E2_V3-E2_V2 -9.000000e+00 -12.8865633 -5.1134367 0.0000350 #E3_V3-E2_V2 -5.500000e+00 -9.3865633 -1.6134367 0.0038367 #E1_V4-E2_V2 -9.000000e+00 -12.8865633 -5.1134367 0.0000350 #E2_V4-E2_V2 -1.200000e+01 -15.8865633 -8.1134367 0.0000016 #E3_V4-E2_V2 -1.100000e+01 -14.8865633 -7.1134367 0.0000041 #E1_V3-E3_V2 -8.000000e+00 -11.8865633 -4.1134367 0.0001171 #E2_V3-E3_V2 -6.000000e+00 -9.8865633 -2.1134367 0.0018056 #E3_V3-E3_V2 -2.500000e+00 -6.3865633 1.3865633 0.3913006 #E1_V4-E3_V2 -6.000000e+00 -9.8865633 -2.1134367 0.0018056 #E2_V4-E3_V2 -9.000000e+00 -12.8865633 -5.1134367 0.0000350 #E3_V4-E3_V2 -8.000000e+00 -11.8865633 -4.1134367 0.0001171 #E2_V3-E1_V3 2.000000e+00 -1.8865633 5.8865633 0.6637770 #E3_V3-E1_V3 5.500000e+00 1.6134367 9.3865633 0.0038367 #E1_V4-E1_V3 2.000000e+00 -1.8865633 5.8865633 0.6637770 #E2_V4-E1_V3 -1.000000e+00 -4.8865633 2.8865633 0.9934341 #E3_V4-E1_V3 -8.881784e-15 -3.8865633 3.8865633 1.0000000 #E3_V3-E2_V3 3.500000e+00 -0.3865633 7.3865633 0.0926317 #E1_V4-E2_V3 -1.776357e-15 -3.8865633 3.8865633 1.0000000 #E2_V4-E2_V3 -3.000000e+00 -6.8865633 0.8865633 0.1985174 #E3_V4-E2_V3 -2.000000e+00 -5.8865633 1.8865633 0.6637770 #E1_V4-E3_V3 -3.500000e+00 -7.3865633 0.3865633 0.0926317 #E2_V4-E3_V3 -6.500000e+00 -10.3865633 -2.6134367 0.0008735 #E3_V4-E3_V3 -5.500000e+00 -9.3865633 -1.6134367 0.0038367 #E2_V4-E1_V4 -3.000000e+00 -6.8865633 0.8865633 0.1985174 #E3_V4-E1_V4 -2.000000e+00 -5.8865633 1.8865633 0.6637770 #E3_V4-E2_V4 1.000000e+00 -2.8865633 4.8865633 0.9934341 # representation graphique : par(mar=c(4, 4, 4, 3)) plot(TukeyHSD(modINTER2)) # approche par specifications de contrastes d interet # # par exemple : #tester Variete V1 contre Variete V2 pour Engrais =E1 # cad E1_V1 contre E1_V2 # et # tester Engrais 1 contre Engrais 3 pour Variete =V2 # cad E1_V2 contre E3_V2 # et # comme E1_V1 n'est pas differente de E1_V3, les reunir en E1_V1+ E1_V3+ et le tester contre E1_V2 #interac #Levels: E1_V1 E2_V1 E3_V1 E1_V2 E2_V2 E3_V2 E1_V3 E2_V3 E3_V3 E1_V4 E2_V4 E3_V4 contraste1 <- c( 1,-1,rep(0,10) ) contraste2 <- c( rep(0,3),1,0,-1,rep(0,6)) contraste3 <- c( 1/2,0,1/2,-1,rep(0,8)) # construire la matrice de contrastes que l on veut tester M<-rbind("E1_V1 contre E1_V2"=contraste1,"E1_V2 contre E3_V2"=contraste2, "E1-V3+E1_V1 contre E_1_V2"=contraste3) M test<-multcomp::glht(modINTER,linfct=mcp(interac=M)) #test<-glht(modINTER2,linfct=mcp(interac=M)) marche aussi et mem resultat ! summary(test) # tous les comparaisons du facteur Varietes pour les differents niveaux d'Engrais mod_mais<-lm(prod~Engrais*Varietes,data=Mais) ############################### # Approche 2 : travailler sur les moyennes ajustées (et approche INDISPENSABLE si expérience non équilibrée) # approche permettant de comparer deux à deux tous les niveaux d'un facteur # conditionnelement aux niveaux de l'autre sur la base des moyennes ajustées library(emmeans) emmeans(mod_mais, specs=pairwise~Varietes|Engrais) # tous les comparaisons du facteur Varietes pour les differents niveaux d'Engrais # avec V1 pris comme ref emmeans(mod_mais, specs=dunnett~Varietes|Engrais) # voir ? emmeans() pour choix possibles de spec=