silico.biotoul.fr
 

M1 Traitement de Donnees Biologiques - TP 4 R

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
 
(46 intermediate revisions not shown)
Line 1: Line 1:
= Tests statistiques:  Analyse de la Variance (ANOVA) =
= Tests statistiques:  Analyse de la Variance (ANOVA) =
-
L'Analyse de la Variance (ANOVA) est une <span style='color: #990000'>extension du test de Student de comparaison de moyennes, à  K>2 groupes (K=2 groupes: test de Student).</span>
+
L'Analyse de la Variance (ANOVA) est une '''extension du test de Student de comparaison de moyennes, à  ''k''>2 groupes (''k''=2 groupes : test de Student).'''
-
Une même variable quantitative est donc mesurée dans  K>2 groupes.
+
Une même variable quantitative est donc mesurée dans  ''k''>2 groupes.
Cependant, avant de se lancer dans une comparaison de groupes 2 à 2 afin de détecter des différences, l'intérêt particulier de l'ANOVA est de tester s'il y a un effet général du (ou des) facteur(s) testé(s).
Cependant, avant de se lancer dans une comparaison de groupes 2 à 2 afin de détecter des différences, l'intérêt particulier de l'ANOVA est de tester s'il y a un effet général du (ou des) facteur(s) testé(s).
-
Les groupes sont donc des modalités (ou niveaux) d'un ou plusieurs facteurs (i.e. variable(s) qualitative(s)).
+
Les groupes sont donc des modalités (ou niveaux) d'un ou plusieurs facteurs (''i.e.'' variable(s) qualitative(s)).
-
L'écart d'un ou plusieurs groupes par rapport à la moyenne générale est estimé par la variance inter-groupe. En tenant compte de la variance intra-groupe, un test de Fischer sur ces 2 estimateurs de variance est effectué.
+
L'écart d'un ou plusieurs groupes par rapport à la moyenne générale est estimé par la variance inter-groupe. En tenant compte de la variance intra-groupe, un test de Fisher sur ces 2 estimateurs de variance est effectué.
-
Le tableau de l'ANOVA donne le résultat de ce test de Fischer pour chaque facteur (variable qualitative) testé.
+
Le tableau de l'ANOVA donne le résultat de ce test de Fisher pour chaque facteur (variable qualitative) testé.
-
Dans ce TP, nous verrons:
+
Dans ce TP, nous verrons :
* l'ANOVA à 1 facteur
* l'ANOVA à 1 facteur
* l'ANOVA à 2 facteurs
* l'ANOVA à 2 facteurs
-
<span style='color: #990000'>ATTENTION: l'ANOVA étant un test statistique, il est là aussi crucial de bien connaître l'hypothèse nulle du test (H0) pour pouvoir correctement interpréter !</span>
+
'''ATTENTION : l'ANOVA étant un test statistique, il est là aussi crucial de bien connaître l'hypothèse nulle du test (H0) pour pouvoir correctement interpréter !'''
-
http://silico.biotoul.fr/site/index.php?title=M1_Traitement_de_Donnees_Biologiques_-_TP_4_R&action=edit
+
-
<span style='color: #990000'>Créez tout d'abord un répertoire de travail sur le bureau (par exemple TDB-TP4) et commencez par télécharger le fichier source que vous allez utiliser et compléter pour générer le compte rendu de TP : [[Media:M1.TDB.TP_tests_ANOVA_R.Rmd|M1.TDB.TP_tests_ANOVA_R.Rmd]] (click droit de la souris -- enregistrer la cible sous...).</span> Ouvrez le logiciel RStudio et chargez ce fichier puis lancez sa compilation pour voir le compte rendu. Pour cela cliquez sur le bouton Knit HTML ou bien utilisez la combinaison de touches Ctrl + shift + K.  
+
<span style='color: #990000'>Créez tout d'abord un répertoire de travail sur le bureau (par exemple TDB-TP4) et commencez par télécharger le fichier source que vous allez utiliser et compléter pour générer le compte rendu de TP : [[Media:M1.TDB.TP_tests_ANOVA_R.Rmd|M1.TDB.TP_tests_ANOVA_R.Rmd]] (click droit de la souris -- enregistrer la cible sous...).</span> Ouvrez le logiciel RStudio et chargez ce fichier puis lancez sa compilation pour voir le compte rendu. Pour cela cliquez sur le bouton Knit HTML ou bien utilisez la combinaison de touches <tt>Ctrl + shift + K</tt>.
-
= ANOVA à un facteur (une variable qualitative) =
+
= ANOVA à 1 facteur (une variable qualitative) =
-
Dans cet exemple, nous allons évaluer si la résistance à un pathogène donné varie pour des plantes de 4 lignées différentes de l'espèce modèle Medicago truncatula.
+
Dans cet exemple, nous allons évaluer si la résistance à un pathogène donné varie pour des plantes de 4 lignées différentes de l'espèce modèle ''Medicago truncatula''.
-
Avant de comparer, les scores moyens de résistance, nous effectuons une ANOVA à 1 facteur pour évaluer s'il y a un effet "lignée" sur le score de résitance.
+
Avant de comparer les scores moyens de résistance, nous effectuons une ANOVA à 1 facteur pour évaluer s'il y a un effet "lignée" sur le score de résistance.
Si c'est le cas, nous comparerons ensuite les lignées 2 à 2 afin de déterminer la(les)quelle(s) se différencie(nt) des autres.  
Si c'est le cas, nous comparerons ensuite les lignées 2 à 2 afin de déterminer la(les)quelle(s) se différencie(nt) des autres.  
-
Téléchargez et sauvegardez le fichier [[Media:resistance.txt|resistance.txt]] dans le répertoire crée pour ce TP (click droit de la souris -- enregistrer la cible sous...).
+
Téléchargez et sauvegardez le fichier [[Media:resistance.txt|resistance.txt]] dans le répertoire créé pour ce TP (click droit de la souris -- enregistrer la cible sous...).
-
== Les données (petits calculs et graphiques)  ==
+
== Exploration du jeu de données ==
-
Lecture
+
Lecture du fichier
<source lang='rsplus'>
<source lang='rsplus'>
-
resistance = read.table("resistance.txt",header=T); resistance
+
resistance = read.table("resistance.txt",header=T,stringsAsFactors=T); head(resistance) ; tail(resistance)
attach(resistance);names(resistance)
attach(resistance);names(resistance)
</source>
</source>
-
Affichage des effectifs des différents groupes (niveaux du facteur):
+
 
 +
Affichage des effectifs des différents groupes (niveaux du facteur) :
<source lang='rsplus'>
<source lang='rsplus'>
table(lignee)
table(lignee)
</source>
</source>
-
Affichage du score moyen par groupe avec la fonction tapply():
+
 
 +
Affichage du score moyen par groupe avec la fonction <tt>tapply()</tt> :
<source lang='rsplus'>
<source lang='rsplus'>
tapply(score,lignee,mean)
tapply(score,lignee,mean)
</source>
</source>
-
... et même un résumé des données par groupe:
+
 
 +
... et même, un résumé des données par groupe :
<source lang='rsplus'>
<source lang='rsplus'>
tapply(score,lignee,summary)
tapply(score,lignee,summary)
</source>
</source>
-
Représentations graphiques des données
+
Représentations graphiques des données : affichage de l'histogramme des scores, tous groupes confondus
-
Affichage de l'histogramme des scores, tout groupes confondus:
+
<source lang='rsplus'>
<source lang='rsplus'>
-
hist(score,breaks=15)
+
hist(score,breaks=15, ylab='effectifs')
-
# un commentaire ?
+
</source>
</source>
-
Affichage de la "boîte à moustache" des scores pour chaque lignée:
+
 
 +
<span style='color: #990000;'>un commentaire ?
 +
</span>
 +
 
 +
Affichage de la distribution des scores pour chaque lignée avec une boîte à moustache :
<source lang='rsplus'>
<source lang='rsplus'>
boxplot(score~lignee,col="green",ylab="score")
boxplot(score~lignee,col="green",ylab="score")
</source>
</source>
-
<span style='color: #990000;'>Ajoutez ces parties à votre compte rendu.</span>
+
<span style='color: #990000;'>Ajoutez ces parties à votre compte rendu.
 +
</span>
== Tests préalables à l'ANOVA ==
== Tests préalables à l'ANOVA ==
Line 60: Line 65:
=== Test d'adéquation à la loi normale (H0 = "les données suivent la loi Normale") ===
=== Test d'adéquation à la loi normale (H0 = "les données suivent la loi Normale") ===
-
<span style='color: #990000;'>Test de Shapiro</span>
+
Test de Shapiro
<source lang='rsplus'>
<source lang='rsplus'>
shapiro.test(score[lignee=="HM008"])
shapiro.test(score[lignee=="HM008"])
Line 66: Line 71:
shapiro.test(score[lignee=="DZA45"])
shapiro.test(score[lignee=="DZA45"])
shapiro.test(score[lignee=="HM013"])
shapiro.test(score[lignee=="HM013"])
-
# remarque: on anticipe un peu, mais on peut aussi tester la normalité des données sur les résidus de l'ANOVA:
 
-
shapiro.test(residuals(aov(score~lignee)))
 
</source>
</source>
-
 
+
<span style='color: #990000;'>Conclusion ?</span>
-
 
+
=== Test d'homogénéité des variances intra-groupes (H0 = "les variances sont égales") ===
=== Test d'homogénéité des variances intra-groupes (H0 = "les variances sont égales") ===
 +
Test de Bartlett de comparaison de plus de 2 variances :
 +
<source lang='rsplus'>
 +
bartlett.test(score~lignee)
 +
</source>
 +
<span style='color: #990000;'>Conclusion ?</span>
-
== ANOVA à 1 facteur ==
+
== ANOVA à 1 facteur (le facteur est la variable "lignee") => H0 :"les moyennes des différents groupes sont égales" ==
-
 
+
On utilise la commande <tt>aov()</tt> :
-
 
+
<source lang='rsplus'>
-
<span style='color: #990000;'>Ajoutez ces parties à votre compte rendu.</span>
+
res = aov(score~lignee)
-
 
+
</source>
-
 
+
-
 
+
-
 
+
-
<span style='color: #990000;'>Test de Shapiro</span>
+
IMPORTANT : le tableau de l'ANOVA s'obtient avec la fonction <tt>summary()</tt> :
<source lang='rsplus'>
<source lang='rsplus'>
 +
summary(res)  # ou summary(aov(score~lignee))
</source>
</source>
 +
<span style='color: #990000;'>
 +
Conclusion ?</span>
-
 
+
On peut visualiser la différence des moyennes de chaque groupe avec la moyenne générale avec :
-
 
+
-
 
+
-
 
+
-
 
+
-
Suite à un croisement effectué en laboratoire, les proportions génotypiques observées (absolues) parmis 100 individus descendants sont:
+
<source lang='rsplus'>
<source lang='rsplus'>
-
freq_obs =c(22,53,25)
+
model.tables(res)
</source>
</source>
-
Question: les proportions observées sont-elles conformes aux proportions mendéliennes? Pour cela il faut effectuer la commande suivante, et répondre aux questions ci-après:
+
Si l'ANOVA détecte un effet significatif du facteur, on peut chercher les inégalités de moyennes : Tests de Student pour comparer les groupes 2 à 2
<source lang='rsplus'>
<source lang='rsplus'>
-
chisq.test(freq_obs,p=freq_exp)
+
t.test(score[lignee=="DZA45"],score[lignee=="HM013"],var.equal=T)
 +
t.test(score[lignee=="DZA45"],score[lignee=="HM008"],var.equal=T)
 +
t.test(score[lignee=="DZA45"],score[lignee=="A17"],var.equal=T)
 +
t.test(score[lignee=="A17"],score[lignee=="HM013"],var.equal=T)
 +
t.test(score[lignee=="HM008"],score[lignee=="HM013"],var.equal=T)
 +
t.test(score[lignee=="A17"],score[lignee=="HM008"],var.equal=T)
</source>
</source>
-
* Quelle est l'hypothèse nulle (H0) du test ? L'hypothèse alternative (H1) ?
+
IMPORTANT : pour limiter le taux de faux positifs lors de '''tests multiples''', et pour faire toutes ces comparaisons avec une seule commande, on fait un test de comparaisons multiples de Student avec correction de la p-valeur :
-
* Qu'indique la p-valeur ?
+
-
* La p-valeur est-elle inférieure ou supérieure au seuil alpha = 5% ?
+
-
 
+
-
Conclusion:
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
== Tests de Chi2 d'indépendance ==
+
-
 
+
-
<span style='color: #990000'>On cherche à évaluer la relation entre deux variables qualitatives à 2 ou > 2 modalités.</span>
+
-
Le test de Chi2 d'indépendance s'effectue sur les effectifs des différentes classes des 2 variables qualitatives.
+
-
 
+
-
Notre exemple: la résistance ou sensibilité à un pathogène est-elle indépendante de l'écotype d'Arabidopsis thaliana ? Autrement dit y-a-t-il indépendance entre phénotype et écotype?
+
-
 
+
-
On va faire un test de Chi2, sur table de contingence des effectifs des différents écotypes pour chaque phénotype qualitatif, comme celle-ci:
+
-
 
+
-
                  Col  Ws Can
+
-
        Résistant 15  5  3
+
-
        Sensible  0  19  16
+
-
 
+
-
Pour cela on crée un tableau de contingence:
+
<source lang='rsplus'>
<source lang='rsplus'>
-
mat=matrix(c(15,5,3,0,19,16),nrow=2,byrow=T, dimnames=list(c("R","S"), c("Col","Ws","Can")))
+
pairwise.t.test(score,lignee,p.adjust.method="bonferroni")
-
mat
+
</source>
</source>
-
Puis on utilise la commande:
+
On peut aussi faire un test de Duncan (l'intérêt étant qu'il propose de classer les groupes en groupes similaires a, b, ''etc.'') :
<source lang='rsplus'>
<source lang='rsplus'>
-
chisq.test(mat, correct =F)
+
# le test nécessite de charger la librairie laercio
 +
library(laercio)
 +
LDuncan(res,conf.level = 0.99)
</source>
</source>
-
Conclusion (posez vous les mêmes questions que dans l'exercice précédent) ?
 
-
 
-
S'il y a une relation entre phénotype et écotype, comment se traduit-elle ?
 
<source lang='rsplus'>
<source lang='rsplus'>
-
obs=chisq.test(mat, correct =F)$observed
+
# Si version R > 4.2, utiliser la fonction duncan.test du package agricolae
-
exp=chisq.test(mat, correct =F)$expected
+
install.packages("agricolae")
-
obs-exp
+
library(agricolae)
 +
duncan.test(res,"lignee",alpha=0.01,console=T)
</source>
</source>
-
Conclusion ?
 
-
<span style='color: #990000;'>Ajoutez ces parties à votre compte rendu.</span>
 
 +
<span style='color: #990000;'>
 +
Conclusion de ces analyses ?</span>
-
= Tests paramétriques d'adéquation et tests d'homogénéité (sur une variable quantitative) =
+
<span style='color: #990000;'>Ajoutez ces parties à votre compte rendu.</span>
-
<span style='color: #990000;'> Nous allons intégrer ces tests dans le contexte d'une même analyse statistique. </span>
+
==<span style='color: #990000;'>REMARQUE IMPORTANTE, si les échantillons ne suivent pas la normalité</span>==
-
 
+
On peut faire une ANOVA non paramétrique avec le test de Kruskall-Wallis :
-
Dans cet exercice, nous allons chercher à comparer les valeurs d'une variable quantitative mesurée dans 2 échantillons 1 et 2, à l'aide d'un test de Student de comparaison de moyennes de 2 échantillons.
+
-
<span style='color: #990000;'>La question posée est: la taille moyenne (à un temps de donné) de plantules d' A. thaliana de génotype "sauvage" (Col_0) et d'un mutant pour le gène X1 (Mut_X1) est-elle la même (hypothèse H0)?</span>
+
-
 
+
-
Rappel de la procédure (vue en cours):
+
-
 
+
-
[[Image:test_homogeneite.jpeg]]
+
-
 
+
-
== Les données (petits calculs et graphiques)  ==
+
<source lang='rsplus'>
<source lang='rsplus'>
-
# Vecteurs des tailles des plantes
+
kruskal.test(score~lignee)
-
Col_0=c(4.30,4.25,3.50,3.35,4.30,3.75,3.55,4.10,3.95,4.55,4.25,3.75,3.85,4.15,3.55,4.75,3.95,3.65)
+
-
Mut_X1=c(3.06,4.05,3.95,3.40,3.80,3.95,3.65,4,3.85,3.95,3.65,3.75,3.4)
+
</source>
</source>
 +
Puis des tests non paramétrique de Wilcoxon/Mann-Whitney, pour comparer les groupes 2 à 2 :
<source lang='rsplus'>
<source lang='rsplus'>
-
# Valeurs moyennes
+
wilcox.test(score[lignee=="HM008"],score[lignee=="A17"])
-
mean(Col_0)
+
wilcox.test(score[lignee=="HM008"],score[lignee=="DZA45"])
-
mean(Mut_X1)
+
</source>
</source>
 +
... et, pour les mêmes raisons que pour le test multiple de Student, faire un test multiple de Wilcoxon/Mann-Whitney :
<source lang='rsplus'>
<source lang='rsplus'>
-
# Représentation graphique des données
+
pairwise.wilcox.test(score,lignee)
-
boxplot(Col_0,Mut_X1,names=c("Col_0","Mut_X1"),col=c("white","darkgreen"),ylab="taille (cm)")
+
pairwise.wilcox.test(jitter(score),lignee,p.adjust.method="bonferroni") # pour résoudre le problème des rangs identiques dans le test, on peut rajouter un bruit aux valeurs
-
abline(h=mean(Col_0),col="black",lty=3,lwd=2)
+
-
abline(h=mean(Mut_X1),col="green",lty=3,lwd=2)
+
-
legend("topright", legend = c("Mut_X1","Col_A" ), text.col = c("darkgreen","black"))
+
</source>
</source>
-
== Test d'adéquation à la loi normale (H0 = "les données suivent la loi Normale") ==
+
= ANOVA à 2 facteurs (deux variables qualitatives) =
-
<span style='color: #990000;'>Test de Shapiro</span>
+
-
<source lang='rsplus'>
+
-
shapiro.test(Col_0)
+
-
shapiro.test(Mut_X1)
+
-
</source>
+
-
Conclusion?
+
Dans cet exemple, nous allons évaluer l'effet du type d'engrais ("A" ou "B") et de la température ambiante ("High" ou "Medium") sur la croissance de plantules en chambre de cultures.
-
 
+
Ici, les 2 facteurs ("engrais", "température") ont 2 niveaux chacun. Il y a donc 4 groupes de plantes à comparer.
-
== Test d'homogénéité des variances (H0 = "les variances sont égales") ==
+
Avant de comparer la croissance moyenne des groupes 2 à 2, nous effectuons une ANOVA à 2 facteurs pour évaluer s'il y a un effet "engrais", "température", des 2, ''etc.'' sur la croissance des plantes.
 +
L'ANOVA est plus intéressante quand on regarde l'effet de plusieurs facteurs.
 +
Téléchargez et sauvegardez le fichier [[Media:croissance.txt|croissance.txt]] dans le répertoire crée pour ce TP (click droit de la souris -- enregistrer la cible sous...).
-
<span style='color: #990000;'>Test de Fischer</span>
+
== Exploration du jeu de données ==
 +
Lecture du fichier
<source lang='rsplus'>
<source lang='rsplus'>
-
var.test(Col_0,Mut_X1)  
+
croissance = read.table("croissance.txt",header=T,stringsAsFactors=T) ; head(croissance)
 +
attach(croissance);names(croissance)
</source>
</source>
-
Conclusion?
+
Représentations graphiques des données, utilisation de <tt>boxplot</tt>.
-
== Test d'homogénéité (de comparaison) des moyennes (H0 = "les moyennes sont égales")==
+
'''Remarque :''' avec la fonction <tt>layout</tt>, on passe une matrice créée à la volée (avec un vecteur qui la rempli par colonnes) de 2 lignes et 2 colonnes qui spécifie à quelle figure (1, 2 ou 3) appartient chaque zone :
 +
{| border=1 cellspacing=0 cellpadding=10
 +
|-
 +
| 1 || 3
 +
|-
 +
| 2 || 3
 +
|}
-
<span style='color: #990000;'>Test de Student bilatéral (H1 = "il y a une différence de moyenne entre Col_0 et Mut_X1")</span>
 
-
<source lang='rsplus'>
 
-
# dans la commande t.test on peut spécifier si les variances des échantillons sont égales (var.equal=T) ou non(var.equal=F)
 
-
# l'hypothèse alternative est "two.sided" par défaut
 
-
t.test(Col_0,Mut_X1,var.equal = T)
 
-
</source>
 
-
<span style='color: #990000;'>Test de Student unilatéral (H1 = "la moyenne de Col_0 > (ou <) à la moyenne de Mut_X1")</span>
 
<source lang='rsplus'>
<source lang='rsplus'>
-
# dans la commande t.test, on va spécifier l'hypothèse alternative: alternative = "greater" ou "less". Quel terme allez-vous mettre?
+
layout( matrix( c(1,2,3,3) , nrow=2, ncol=2, byrow = FALSE))
-
t.test(Col_0,Mut_X1,var.equal = T, alternative = "???")
+
boxplot(hauteur~engrais, main="engrais",ylab="hauteur (cm)")
 +
boxplot(hauteur~temperature, main="temperature",ylab="hauteur (cm)")
 +
boxplot(hauteur~temperature+engrais, main="engrais-temperature",ylab="hauteur (cm)",las=3)
</source>
</source>
-
 
-
Conclusion finale?
 
<span style='color: #990000;'>Ajoutez ces parties à votre compte rendu.</span>
<span style='color: #990000;'>Ajoutez ces parties à votre compte rendu.</span>
 +
== Tests préalables à l'ANOVA ==
 +
=== Test d'adéquation à la loi normale (H0 = "les données suivent la loi Normale") ===
 +
Test de Shapiro
 +
<source lang='rsplus'>
 +
# remarque: pour sélectionner les lignes correspondant à une combinaison de facteur colonnes, on utilise "&":
 +
shapiro.test(hauteur[engrais=="A" & temperature=="Medium"])
 +
shapiro.test(hauteur[engrais=="A" & temperature=="High"])
 +
shapiro.test(hauteur[engrais=="B" & temperature=="Medium"])
 +
shapiro.test(hauteur[engrais=="B" & temperature=="High"])
 +
</source>
-
= Test non paramétriques d'homogénéité (sur une variable quantitative) =
+
<span style='color: #990000;'>
 +
Conclusion ?
 +
</span>
-
Si les données ne suivent pas la loi Normale dans chaque échantillon, nous pouvons effectuer le test non paramétrique de Mann & Whitney
+
=== Test d'homogénéité des variances intra-groupes (H0 = "les variances sont égales") ===
 +
Test de Bartlett de comparaison de plus de 2 variances
 +
<source lang='rsplus'>
 +
bartlett.test(hauteur~interaction(temperature,engrais))
 +
</source>
-
<span style='color: #990000;'> Dans cet exercice, nous cherchons à comparer l'expression d'un gène de défense chez la plante en condition d'infection par un microorganisme pathogène foliaire, avec la condition "contrôle" (plante non infectée).
+
<span style='color: #990000;'>Conclusion ?
 +
</span>
-
== Les données ==
+
== ANOVA à 2 facteurs (les facteurs sont "engrais" et "température") => H0 : "les moyennes des différents groupes sont égales" ==
-
Les données sont les valeurs d'expression normalisée du gène pour 3 réplicats techniques sur 3 réplicats biologiques de tissus foliaire.
+
=== Modèle d'ANOVA où l'on teste l'effet de 2 facteurs ainsi que leur interaction ===
 +
<source lang='rsplus'>
 +
model1<-aov(hauteur~temperature*engrais)
 +
summary(model1)
 +
interaction.plot(temperature,engrais,hauteur, main="interaction plot")
 +
</source>
 +
=== Modèle d'ANOVA où l'on teste l'effet de 2 facteurs sans leur interaction ===
<source lang='rsplus'>
<source lang='rsplus'>
-
non_infected=c(0.021,0.15,0.023,0.03,0.022,0.05,0.035,0.1,0.03)
+
model2<-aov(hauteur~temperature+engrais)
-
infected=c(1.22,1.12,1.06,1.04,0.86,1.24,1.96,0.9,2.5)
+
summary(model2)
-
mean(non_infected)
+
-
mean(infected)
+
</source>
</source>
-
== Test d'adéquation à la loi normale  ==
+
=== Comparaison des 2 modèles avec une ANOVA sur les résidus des 2 modèles ===
<source lang='rsplus'>
<source lang='rsplus'>
-
# Test de Shapiro
+
anova(model1,model2)
-
shapiro.test(infected)
+
-
shapiro.test(non_infected)
+
</source>
</source>
-
== Test non paramétrique d'homogénéité ==
+
* s'il n'y a pas de différence (p-value > 0.05) on garde le modèle le plus simple.
 +
* s'il y a une différence significative (p-value < 0.05), on garde le modèle présentant le plus de termes significatifs.
 +
 
 +
<span style='color: #990000;'>
 +
Conclusion?
 +
</span>
 +
 
 +
On peut maintenant comparer les groupes de plantes 2 à 2 :
<source lang='rsplus'>
<source lang='rsplus'>
-
# Test de Mann & Whitney (équivalent non paramétrique du test de Student)
+
# plantes avec "engrais" différents
-
wilcox.test(infected,non_infected)                       # bilatéral
+
t.test(hauteur[engrais=="A" & temperature=="High"],hauteur[engrais=="B" & temperature=="High"],var.equal=T)     # comparaison A/B pour température "High"
-
wilcox.test(infected,non_infected,alternative="greater") # unilatéral
+
t.test(hauteur[engrais=="A" & temperature=="Medium"],hauteur[engrais=="B" & temperature=="Medium"],var.equal=T) # comparaison A/B pour température "Medium"
 +
# plantes avec "températures" différentes
 +
t.test(hauteur[engrais=="A" & temperature=="Medium"],hauteur[engrais=="A" & temperature=="High"],var.equal=T)  # comparaison High/Medium pour l'engrais A
 +
t.test(hauteur[engrais=="B" & temperature=="Medium"],hauteur[engrais=="B" & temperature=="High"],var.equal=T)   # comparaison High/Medium pour l'engrais B
</source>
</source>
 +
<span style='color: #990000;'>
Conclusion?
Conclusion?
 +
</span>
 +
 +
<span style='color: #990000;'>Ajoutez ces parties à votre compte rendu.</span>
 +
 +
= PROBLEME SUPPLEMENTAIRE (à faire à la maison) =
 +
 +
Comparez la production journalière de lait pour 3 races de vaches laitières. Les données sont dans le fichier suivant: [[Media:vache.txt|vache.txt]]
 +
-
<span style='color: #990000;'>Ajoutez tout cela au compte rendu de TP avant de l'envoyer à votre enseignant par mail ([mailto:bonhomme@lrsv.ups-tlse.fr bonhomme@lrsv.ups-tlse.fr] ou [mailto:barriot@biotoul.fr barriot@biotoul.fr]). Le compte rendu est à envoyer '''avant de commencer le TP5'''. Envoyez les 2 fichiers (.Rmd et .html). Envoyez-vous aussi le mail en copie pour pouvoir vérifier que tout est bien passé. Mettez un titre tel que "Compte rendu TP4 TDB de -et votre Nom et Prénom-".</span>
+
<span style='color: #990000;'>Ajoutez tout cela au compte rendu de TP avant de l'envoyer à votre enseignant par mail ([mailto:maxime.bonhomme@univ-tlse3.fr maxime.bonhomme@univ-tlse3.fr] ou [mailto:roland.barriot@univ-tlse3.fr roland.barriot@univ-tlse3.fr]). Le compte rendu est à envoyer '''avant de commencer le TP5'''. Envoyez les 2 fichiers (.Rmd et .html). Envoyez-vous aussi le mail en copie pour pouvoir vérifier que tout est bien passé. Mettez un titre tel que "Compte rendu TP4 TDB de -et votre Nom et Prénom-". '''Votre compte-rendu ne doit pas être un simple "copié-collé" des commandes. Nous attendons une mise en forme correcte (titres, sous-titres, paragraphes, ...), des commentaires (notamment sur les résultats des commandes), des essais de votre part pour améliorer par exemple la visualisation d'un graphique, ou pour rechercher des fonctions R alternatives ou complémentaires à ce qui est proposé; bref une personnalisation de votre compte-rendu !'''</span>
= Liens =
= Liens =

Current revision as of 06:45, 22 September 2022

Contents

Tests statistiques: Analyse de la Variance (ANOVA)

L'Analyse de la Variance (ANOVA) est une extension du test de Student de comparaison de moyennes, à k>2 groupes (k=2 groupes : test de Student). Une même variable quantitative est donc mesurée dans k>2 groupes. Cependant, avant de se lancer dans une comparaison de groupes 2 à 2 afin de détecter des différences, l'intérêt particulier de l'ANOVA est de tester s'il y a un effet général du (ou des) facteur(s) testé(s). Les groupes sont donc des modalités (ou niveaux) d'un ou plusieurs facteurs (i.e. variable(s) qualitative(s)). L'écart d'un ou plusieurs groupes par rapport à la moyenne générale est estimé par la variance inter-groupe. En tenant compte de la variance intra-groupe, un test de Fisher sur ces 2 estimateurs de variance est effectué. Le tableau de l'ANOVA donne le résultat de ce test de Fisher pour chaque facteur (variable qualitative) testé.

Dans ce TP, nous verrons :

  • l'ANOVA à 1 facteur
  • l'ANOVA à 2 facteurs

ATTENTION : l'ANOVA étant un test statistique, il est là aussi crucial de bien connaître l'hypothèse nulle du test (H0) pour pouvoir correctement interpréter !

Créez tout d'abord un répertoire de travail sur le bureau (par exemple TDB-TP4) et commencez par télécharger le fichier source que vous allez utiliser et compléter pour générer le compte rendu de TP : M1.TDB.TP_tests_ANOVA_R.Rmd (click droit de la souris -- enregistrer la cible sous...). Ouvrez le logiciel RStudio et chargez ce fichier puis lancez sa compilation pour voir le compte rendu. Pour cela cliquez sur le bouton Knit HTML ou bien utilisez la combinaison de touches Ctrl + shift + K.

ANOVA à 1 facteur (une variable qualitative)

Dans cet exemple, nous allons évaluer si la résistance à un pathogène donné varie pour des plantes de 4 lignées différentes de l'espèce modèle Medicago truncatula. Avant de comparer les scores moyens de résistance, nous effectuons une ANOVA à 1 facteur pour évaluer s'il y a un effet "lignée" sur le score de résistance. Si c'est le cas, nous comparerons ensuite les lignées 2 à 2 afin de déterminer la(les)quelle(s) se différencie(nt) des autres. Téléchargez et sauvegardez le fichier resistance.txt dans le répertoire créé pour ce TP (click droit de la souris -- enregistrer la cible sous...).

Exploration du jeu de données

Lecture du fichier

resistance = read.table("resistance.txt",header=T,stringsAsFactors=T); head(resistance) ; tail(resistance)
attach(resistance);names(resistance)

Affichage des effectifs des différents groupes (niveaux du facteur) :

table(lignee)

Affichage du score moyen par groupe avec la fonction tapply() :

tapply(score,lignee,mean)

... et même, un résumé des données par groupe :

tapply(score,lignee,summary)

Représentations graphiques des données : affichage de l'histogramme des scores, tous groupes confondus

hist(score,breaks=15, ylab='effectifs')

un commentaire ?

Affichage de la distribution des scores pour chaque lignée avec une boîte à moustache :

boxplot(score~lignee,col="green",ylab="score")

Ajoutez ces parties à votre compte rendu.

Tests préalables à l'ANOVA

Avant de procéder à l'ANOVA, il faut vérifier la normalité des données de chaque groupe, et l'homogénéité des variances des groupes:

Test d'adéquation à la loi normale (H0 = "les données suivent la loi Normale")

Test de Shapiro

shapiro.test(score[lignee=="HM008"])
shapiro.test(score[lignee=="A17"] )
shapiro.test(score[lignee=="DZA45"])
shapiro.test(score[lignee=="HM013"])

Conclusion ?

Test d'homogénéité des variances intra-groupes (H0 = "les variances sont égales")

Test de Bartlett de comparaison de plus de 2 variances :

bartlett.test(score~lignee)

Conclusion ?

ANOVA à 1 facteur (le facteur est la variable "lignee") => H0 :"les moyennes des différents groupes sont égales"

On utilise la commande aov() :

res = aov(score~lignee)

IMPORTANT : le tableau de l'ANOVA s'obtient avec la fonction summary() :

summary(res)  # ou summary(aov(score~lignee))

Conclusion ?

On peut visualiser la différence des moyennes de chaque groupe avec la moyenne générale avec :

model.tables(res)

Si l'ANOVA détecte un effet significatif du facteur, on peut chercher les inégalités de moyennes : Tests de Student pour comparer les groupes 2 à 2

t.test(score[lignee=="DZA45"],score[lignee=="HM013"],var.equal=T)
t.test(score[lignee=="DZA45"],score[lignee=="HM008"],var.equal=T)
t.test(score[lignee=="DZA45"],score[lignee=="A17"],var.equal=T)
t.test(score[lignee=="A17"],score[lignee=="HM013"],var.equal=T)
t.test(score[lignee=="HM008"],score[lignee=="HM013"],var.equal=T)
t.test(score[lignee=="A17"],score[lignee=="HM008"],var.equal=T)

IMPORTANT : pour limiter le taux de faux positifs lors de tests multiples, et pour faire toutes ces comparaisons avec une seule commande, on fait un test de comparaisons multiples de Student avec correction de la p-valeur :

pairwise.t.test(score,lignee,p.adjust.method="bonferroni")

On peut aussi faire un test de Duncan (l'intérêt étant qu'il propose de classer les groupes en groupes similaires a, b, etc.) :

# le test nécessite de charger la librairie laercio
library(laercio)
LDuncan(res,conf.level = 0.99)
# Si version R > 4.2, utiliser la fonction duncan.test du package agricolae
install.packages("agricolae")
library(agricolae)
duncan.test(res,"lignee",alpha=0.01,console=T)


Conclusion de ces analyses ?

Ajoutez ces parties à votre compte rendu.

REMARQUE IMPORTANTE, si les échantillons ne suivent pas la normalité

On peut faire une ANOVA non paramétrique avec le test de Kruskall-Wallis :

kruskal.test(score~lignee)

Puis des tests non paramétrique de Wilcoxon/Mann-Whitney, pour comparer les groupes 2 à 2 :

wilcox.test(score[lignee=="HM008"],score[lignee=="A17"])
wilcox.test(score[lignee=="HM008"],score[lignee=="DZA45"])

... et, pour les mêmes raisons que pour le test multiple de Student, faire un test multiple de Wilcoxon/Mann-Whitney :

pairwise.wilcox.test(score,lignee)
pairwise.wilcox.test(jitter(score),lignee,p.adjust.method="bonferroni") # pour résoudre le problème des rangs identiques dans le test, on peut rajouter un bruit aux valeurs

ANOVA à 2 facteurs (deux variables qualitatives)

Dans cet exemple, nous allons évaluer l'effet du type d'engrais ("A" ou "B") et de la température ambiante ("High" ou "Medium") sur la croissance de plantules en chambre de cultures. Ici, les 2 facteurs ("engrais", "température") ont 2 niveaux chacun. Il y a donc 4 groupes de plantes à comparer. Avant de comparer la croissance moyenne des groupes 2 à 2, nous effectuons une ANOVA à 2 facteurs pour évaluer s'il y a un effet "engrais", "température", des 2, etc. sur la croissance des plantes. L'ANOVA est plus intéressante quand on regarde l'effet de plusieurs facteurs. Téléchargez et sauvegardez le fichier croissance.txt dans le répertoire crée pour ce TP (click droit de la souris -- enregistrer la cible sous...).

Exploration du jeu de données

Lecture du fichier

croissance = read.table("croissance.txt",header=T,stringsAsFactors=T) ; head(croissance)
attach(croissance);names(croissance)

Représentations graphiques des données, utilisation de boxplot.

Remarque : avec la fonction layout, on passe une matrice créée à la volée (avec un vecteur qui la rempli par colonnes) de 2 lignes et 2 colonnes qui spécifie à quelle figure (1, 2 ou 3) appartient chaque zone :

1 3
2 3


layout( matrix( c(1,2,3,3) , nrow=2, ncol=2, byrow = FALSE))
boxplot(hauteur~engrais, main="engrais",ylab="hauteur (cm)")
boxplot(hauteur~temperature, main="temperature",ylab="hauteur (cm)")
boxplot(hauteur~temperature+engrais, main="engrais-temperature",ylab="hauteur (cm)",las=3)

Ajoutez ces parties à votre compte rendu.

Tests préalables à l'ANOVA

Test d'adéquation à la loi normale (H0 = "les données suivent la loi Normale")

Test de Shapiro

# remarque: pour sélectionner les lignes correspondant à une combinaison de facteur colonnes, on utilise "&":
shapiro.test(hauteur[engrais=="A" & temperature=="Medium"])
shapiro.test(hauteur[engrais=="A" & temperature=="High"])
shapiro.test(hauteur[engrais=="B" & temperature=="Medium"])
shapiro.test(hauteur[engrais=="B" & temperature=="High"])

Conclusion ?

Test d'homogénéité des variances intra-groupes (H0 = "les variances sont égales")

Test de Bartlett de comparaison de plus de 2 variances

bartlett.test(hauteur~interaction(temperature,engrais))

Conclusion ?

ANOVA à 2 facteurs (les facteurs sont "engrais" et "température") => H0 : "les moyennes des différents groupes sont égales"

Modèle d'ANOVA où l'on teste l'effet de 2 facteurs ainsi que leur interaction

model1<-aov(hauteur~temperature*engrais)	
summary(model1)
interaction.plot(temperature,engrais,hauteur, main="interaction plot")

Modèle d'ANOVA où l'on teste l'effet de 2 facteurs sans leur interaction

model2<-aov(hauteur~temperature+engrais)	
summary(model2)

Comparaison des 2 modèles avec une ANOVA sur les résidus des 2 modèles

anova(model1,model2)
  • s'il n'y a pas de différence (p-value > 0.05) on garde le modèle le plus simple.
  • s'il y a une différence significative (p-value < 0.05), on garde le modèle présentant le plus de termes significatifs.

Conclusion?

On peut maintenant comparer les groupes de plantes 2 à 2 :

# plantes avec "engrais" différents
t.test(hauteur[engrais=="A" & temperature=="High"],hauteur[engrais=="B" & temperature=="High"],var.equal=T)     # comparaison A/B pour température "High"
t.test(hauteur[engrais=="A" & temperature=="Medium"],hauteur[engrais=="B" & temperature=="Medium"],var.equal=T) # comparaison A/B pour température "Medium"
# plantes avec "températures" différentes
t.test(hauteur[engrais=="A" & temperature=="Medium"],hauteur[engrais=="A" & temperature=="High"],var.equal=T)   # comparaison High/Medium pour l'engrais A
t.test(hauteur[engrais=="B" & temperature=="Medium"],hauteur[engrais=="B" & temperature=="High"],var.equal=T)   # comparaison High/Medium pour l'engrais B

Conclusion?

Ajoutez ces parties à votre compte rendu.

PROBLEME SUPPLEMENTAIRE (à faire à la maison)

Comparez la production journalière de lait pour 3 races de vaches laitières. Les données sont dans le fichier suivant: vache.txt



Ajoutez tout cela au compte rendu de TP avant de l'envoyer à votre enseignant par mail (maxime.bonhomme@univ-tlse3.fr ou roland.barriot@univ-tlse3.fr). Le compte rendu est à envoyer avant de commencer le TP5. Envoyez les 2 fichiers (.Rmd et .html). Envoyez-vous aussi le mail en copie pour pouvoir vérifier que tout est bien passé. Mettez un titre tel que "Compte rendu TP4 TDB de -et votre Nom et Prénom-". Votre compte-rendu ne doit pas être un simple "copié-collé" des commandes. Nous attendons une mise en forme correcte (titres, sous-titres, paragraphes, ...), des commentaires (notamment sur les résultats des commandes), des essais de votre part pour améliorer par exemple la visualisation d'un graphique, ou pour rechercher des fonctions R alternatives ou complémentaires à ce qui est proposé; bref une personnalisation de votre compte-rendu !

Liens