ANOVA

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

0 - Téléchargez libraries et les données dans R studio

library(MASS)
library(ggplot2)
library(lattice)

Charger les données dans l’environnement de travail.

data("birthwt")
dim(birthwt)

Les données ont été recueillies au Baystate Medical Center, Springfield, Mass, en 1986. La base de données birthwt compte 189 lignes et 10 colonnes.

1 - Explorez les données

Observez les données. La commande head() permet d’afficher les premières observations d’une variable.

head(birthwt,10)

Les données sont multivariés car représentées par des variables organisées en colonnes et des observations figurant en ligne. Chaque ligne du fichier représente donc une unité statistique sur laquelle on a collecté plusieurs mesures ou données. Ici il y a 10 colonnes soit dix variables.

Observez les noms des variables. Les dix colonnes correspondent aux variables suivantes:

variable.names(birthwt)

Les variables sont: 1. low : statut pondéral 2. age : age de la mère 3. lwt : poids de la mère 4. race : catégorie de couleur de peau 5. smoke : consommation de tabac 6. plt : nombre d’accouchements anteérieurs 7. ht : hypertension 8. ui : irritabilité utérine 9. ftv : consultations gynécologiques 10. bwt : poids des bébés à la naissance

Renommez une variable pour faciliter l’interprétation des données.

Cette variable est l’ethnicité de la mère qui était codé en trois classes 1,2,3 et qu’on peut remplacer par “White”,“Black”,“Other”.

birthwt$race <- factor(birthwt$race,levels=c(1,2,3), 
                       labels=c("White","Black","Other"))
birthwt$race 

On peut regarder les cinq premières valeurs pour le poids des bébés.

birthwt$bwt[1:5]

Certaines variables sont binaires comme low, smoke, ht et ui , tandis que certains variables sont numériques avec soit des valeurs discrètes comme ftv, soit des valeurs continues comme lwt et bwt.

On pourrait se demander quel est le poids moyen des femmes qui fumait durant leur grossesse.

mean(birthwt$lwt[birthwt$smoke==1])

On peut sauvegarder ces données dans un fichier externe.

write.table(birthwt, "birthwt.txt", sep="/t", quote=FALSE, row.names=FALSE)

2 - Les hypothèses derrière l’ANOVA

Observez la distribution des variables.

summary(birthwt$bwt)

Explorez la distribution des données avec des fonctions simples.

boxplot(birthwt)

boxplot(birthwt$lwt, xlab="Poids des bebes")

hist(birthwt$lwt, xlab="Poids des bebes")

Condition 1 de l’ANOVA : Testez la normalité des données sur la colonne bwt.

shapiro.test(birthwt$bwt)

Dans l’exemple ci-dessus, la p-value est non significative (p-value = 0.4353). L’échantillon suit donc une loi normale.

Quelle que soit l’option, tous les tests de Shapiro qui présentent une pvalue > 0.05, donc au risque de alpha = 5%, indiquent qu’on ne peut pas rejeter H0. H0 étant : les données ne suivent pas la loi normale. H1 étant : les données suivent la loi normale. Donc le poids des bébés suit bien une loi Normale: OK pour la première condition.

Condition 2 de l’ANOVA : Testez l’égalité des variance entre race et bwt.

bartlett.test(bwt ~ smoke, birthwt)

La pvalue obtenue est de 0.2196 donc est > 0.05, donc, au risque de alpha = 5%, je ne peux pas rejeter H0, donc le poids des bébés et la race ont des variances qui ne sont pas significativement différentes, donc elles sont égales : OK pour la deuxième condition.

Hypothèse 3 de l’ANOVA : independance des données. Le poids des bébés constituant les groupes de race doivent être indépendants. OK pour la troisième condition car c’est en accord avec le protocole d’échantillonnage (chaque poids correspond à un bébé unique).

2 - L’ANOVA à un facteur

Est ce que le poids des bébés à la naissance diffère-t-il selon si la mère est fumeuse ou non?

Le tableau de données birthwt contient tout pour faire l’anova avec la fonction aov()è. La formule s’écrit de maniere générale : Variable_reponse ~ Facteur1

On va stocker le résultat du modele de l’ANOVA dans une variable bwt_smoke qu’on peut ainsi réutiliser.

bwt_smoke <- aov(bwt ~ smoke, birthwt)
bwt_smoke

On voit que le tableau contient uniquement une estimation de la variance résiduelle (Residual Standar Error) et la somme des carrés (Deg. of Freedom).

Pour afficher le tableau d’analyse de variance il faut utiliser la commande summary().

summary(bwt_smoke)

Une réprésentation graphique pertinente pour ce type de données (une variable numérique décrite par une variable quantitative) est un diagramme de type boîte à moustaches, pour avoir une idée de la variation à l’intérieur des groupes.

bwplot(race~bwt, data=birthwt)

On pourrait aussi utiliser ggplotpour représenter nos données à l’aide de la fonction geom_boxplot. Cette fonction nous permet de représenter des boîtes à moustache en indiquant en y la variable numérique, et en x la variable qualitative contenant les classes qu’on souhaite comparer. Ainsi, si on veut comparer la répartition du poids des bébés en fonction de l’ethnicité de la mère, on pourra faire la commande suivante.

ggplot(data=birthwt, aes(y=bwt, x=factor(smoke)))+
  geom_boxplot()+
  theme_classic()

Observer les moyennes par groupe. La commande aggregate renvoit ses résultats sous la forme d’un dataframe ce qui parfois peut servir pour des analyses ultérieures.

L’usage de cette commande est facile : il suffit d’indiquer le type de relation à laquelle on s’interesse entre deux variables, ici le poids des bébés en fonction de l’ethnicité de la mère, le jeu de données utilisé, ici le birthwt, et la statistique qu’on souhaite obtenir, ici la moyenne.

bwtmeans <- aggregate(bwt~ race, data=birthwt, mean)
bwtmeans

3 - Les comparaisons multiples

Pour comparer plus d’une paires de moyennes il existe des techniques de comparaisons multiples de type Bonferroni. La correction de Bonferroni est une méthode utilisée en statistique pour corriger le seuil de significativité lors de comparaisons multiples.

Ici on va utiliser la variable race car il y a trois groupes with(birthwt, pairwise.t.test(bwt,race,p.adjust.method = “bonferroni”)).

La commande pairwise.test permet de tester toutes les paires moyennes à l’aide du test de Student.

t1 <- t.test(bwt ~ race, data=birthwt,subset=race!="Black")
t2 <- t.test(bwt ~ race, data=birthwt,subset=race!="Other")

Multipliez les P-valeurs obtenue depuis le test de student avec nombre de comparaisons que vous effectuez. Ici on a trois groupes donc trois comparaisons.

c(t1$p.value, t2$p.value)*3

Après multiplications de ces P-valeurs, sont-elles bien égales à la P-valeur obtenue avec la correction de Bonferroni ?

with(birthwt, pairwise.t.test(bwt,race,p.adjust.method = "bonferroni",
                              pool.sd=FALSE))

Pour comparer les deux approches, il est nécessaire de préciser à R pool.sd=FALSE pour utilise la variance propres aux seuls groupes pris en compte dans la comparaison (et non l’ensemble de la variance).

4 - L’ANOVA à 2 facteurs

La formule s’écrit de manière générale : `Variable_réponse ~ Facteur1 + Facteur2

Est ce le poids des nouveaux bébés est influencé par le poids de la mère et son ehnicité

Tester un modèle d’ANOVA à deux facteurs.

bwt_race_ht <- aov(bwt ~ race + ht, birthwt)

Construire un tableau d’ANOVA à deux facteurs.

summary(aov(bwt ~ race + ht, birthwt))

Tester un modèle d’ANOVA à deux facteurs avec intéraction. En effet, dans un modèle statistique, il y a l’hypothèse implicite que chaque variable explicative est indépendante des autres. Cependant, cela n’est pas toujours le cas et par exemple, l’effet de la couleur de peau peut varier sur l’hypertension artérielle. Il est dès lors nécessaire de prendre en compte dans son modèle les effets d’interaction entre les variables.

bwt_race_ht <- aov(bwt ~ race + ht + race:ht, birthwt)

Construire un tableau d’ANOVA d’un modèle avec interaction.

summary(bwt_race_ht)

Obtenir les moyennes d’un modèle d’ANOVA à deux facteurs.

aggregate(bwt ~ race+ht, data=birthwt, mean)

5 - La corrélation

*** Est ce que le poids des nouveaux nés est corrélé a celui des mères?

Statistique descriptive sur le poids des nouveaux nés (bwt) et le poids des mères (lwt).

summary(birthwt[,c("bwt","lwt")])

Visualiser directement les histogrammes de ces deux variables.

histogram(~bwt+lwt, data=birthwt, breaks=NULL, xlab="Histogrammes des poids des meres
          \n et des nouveaux-nes", 
          scales=(x=list(relation="free")))

On peut aussi visualiser les variables avec des courbes de densités en utilisant la fonction densityplot.

densityplot(~age | low, data=birthwt, type="count", bwt=2.5)

Observer la relation entre ces deux variables avant d’effectuer un test statistique. *** Utiliserez vous un test paramétrique ou non paramétrique?***

birthwt$lwt <- birthwt$lwt/2.2
xyplot(bwt~lwt, data=birthwt)

Test paramétrique: mesurer le coefficient de corrélation de Pearson.

with(birthwt, cor(bwt, lwt))

Mesurer le coefficient de corrélation en considerant les paires complètes en excluant des données manquantes.

with(birthwt, cor(bwt, lwt, method="spearman"))

Calculer cette corrélation sur plus de deux variables.

cor(birthwt[c("bwt", "lwt","age")])

Quelle paire de variable montre la plus importante corrélation

Calculer les intervalles de confiance avec la fonction cor.test.

with(birthwt, cor(bwt, lwt))

Restreindre le calcul des intervalles de confiance avec l’option subset, ici on indique qu’on souhaite considérer uniquement les bébés avec un poids au dessus de 2500 grammes.

cor.test(~ bwt+lwt, data=birthwt, subset=bwt>2500)

Test non-paramétrique: mesurer le coefficient de corrélation de Spearman. Ajouter l’option method=“spearman” pour realiser un test d’inference sur les rangs. Le coefficient de Spearman permet de détecter des tendances monotones (fonction monotone soit croissante, soit décroissante). Lorsque la tendance est linéaire, il se comporte de façon similaire au coefficient de Pearson. En revanche, il sera plus élevé que la corrélation de Pearson si la tendance est monotone. Plus la tendance monotone est marquée, plus la valeur du coefficient est proche de 1 ou -1.

cor.test(~ bwt+lwt, data=birthwt, subset=bwt>2500, method="spearman")

6 - La régression linéaire

Modéliser la relation entre le poids des nouveaux nés (variable réponse) et le poids des mères (variable explicative).

lm(bwt ~ lwt, data=birthwt)

Pour obtenir le tableau des coefficients de regression, et leur erreur standard et leur degré de significativité.

bwt_lwt <- lm(bwt ~ lwt, data=birthwt)
summary(bwt_lwt)

Visualiser la droite de régression sur le diagramme de dispersion.

xyplot(bwt ~ lwt, data=birthwt, type=c("p","g","r"))

Tester la significativité du modèle.

anova(bwt_lwt)

Prédiction a partir du modèle. On cherche a calculer les valeurs prédites des poids des bébés en fonction du poids de la mère.

head(fitted(bwt_lwt))

Comparer les valeurs prédites aux valeurs observées.

head(cbind(birthwt$bwt, fitted(bwt_lwt)))

Ajouter au diagramme de dispersion les valeurs prédites.

xyplot(bwt + fitted(bwt_lwt) ~ lwt, data=birthwt)

Estimer les valeurs prédites par le modèle.

dp <- data.frame(lwt=seq(35,120,by=5))
bwtpred <- predict(bwt_lwt, newdata = dp)
dp$bwt <- bwtpred
head(dp)

7 - L’analyse des residus

Utiliser la commande resid() pour obtenir les valeurs des résidus du modèle de regression.

head(resid(bwt_lwt))

Visualiser la dispersion des résidus.

xyplot(resid(bwt_lwt) ~ fitted(bwt_lwt), type=c("p","g"))

8 - Exercices

Exercice 1:

Un chercheur souhaite évaluer l’effet de la concentration (mg/mL) en un nutriment X sur la croissance de bacteries (DO). La croissance est mesurée par densité optique.

Ouvrir les fichiers contenant les données.

exo2<-read.table(file = "bacteries.txt",sep=" ", header = T)

Explorez les fichiers.

dim(exo2)
hist(exo2$concentration)

hist(exo2$DO)

plot(exo2$concentration,exo2$DO)

Il y a n=20 pour chaque variable quantitatives continues, mais les histogrammes ne presentent pas une courbe en forme de cloche comme pour la distribution normale.

Tenter de faire une regression linéaire pour expliquer la variation de croissance bactérienne par la variation de la concentration en nutriment. Vérifier la condition d’application de la régression.

Hypothese 1 : Test de normalite.

shapiro.test(exo2$concentration)
shapiro.test(exo2$DO)

Dans les deux tests avec le risque de 5%, la pvalue > 0.05 donc je ne peux pas rejeter H0, les donnees suivent une loi Normale, je peux donc tenter une regression de type y=ax+b.

modele<-lm(exo2$DO~exo2$concentration)
summary(modele)

L’ordonnée a l’origine (b) et la pente (a) présentent toutes les deux des valeurs estimées (b=104.0633 et a=10.1154) qui sont significativement differentes de 0.

Avant de garder ce modèle je dois vérifier pour les residus leur normalité et leur homoscedasticité

shapiro.test(modele$residuals)
plot(modele$residuals~exo2$concentration)

Que se passe-t-il?

Les résidus ne suivent pas une loi normale donc on devrait s’arreter la et tenter soit de transformer l’une des variables, soit de prendre un autre modèle que le modèle linéaire simple, soit chercher à obtenir plus de données de concentration bactérienne en fonction des nutriments. Regardons quand meme visuellement les résidus avec la fonction plot interactive.

plot(modele)

Le 2eme graphe nous montre que les résidus ne sont pas bien alignés sur la droite de normalité. Le 4eme graphe nous montre que le point #1 est au dela de l’isoligne 0.5. C’est sans doute ce point qui peut induire la non normalité des residus car seulement nous n’avons que 20 données, donc cet ecart à la droite de régression peut avoir beaucoup d’impact. Les autres graphes sont satisfaisants.

Avec ce genre d’analyse des résidus, et vu que le coefficient de Determination R^2 vaut 0.9143, donc que la part de variance de Y expliquée par la variation de x est de 91% - ce qui est très bien - on peut presque dire que notre modèle est bon. On peut également prédire qu’il sera possibel qu’il le devienne vraiment statistiquement parlant avec quelques mesures de plus qui rendront certainement les residus normalement distribués.

Pour confirmer cette prédiction, on peut ajouter la droite du modèle trouvé.

plot(exo2$concentration,exo2$DO)
abline(a=modele$coefficients[[1]],b=modele$coefficients[[2]], col="red")

Exercice 2:

Une équipe d’ornithologues a suivi pendant tout un été des nichoirs et a systematiquement dénombré et mesuré la longueur du tarse (pattes) et le poids des poussins.

Quelle question biologique se sont-ils posés?

Reponse : le poids et le tarse des poussins covarient ils? c’est a dire est ce que le poids et la longueur du tars sont corrélés?

Ici on ne cherche pas expliquer la variation de l’un par la variation de l’autre, on veut juste regarder s’ils varient ou non dans le même sens ou le sens opposé. L’hypothèse biologique serait une covariation positive. On va donc faire une analyse de corrélation.

J’ouvre le fichier, ici les colonnes sont separées par des tabulations donc on le specifie avec l’option sep.

exo3=read.table("mesures_nichoirs.txt", sep="\t", header=TRUE)
head(exo3)

Observer la structure du tableau de données.

str(exo3)

Transformer la variable nichoir en variable categorielle.

exo3$nichoir.fac<-factor(exo3$nichoir)

Sous echantilloner le jeu de données avec seulement l’espèce PARCAE.

CAE<-subset(exo3, espece=="PARCAE")

Explorer la distribution des variables.

plot(tarse~poids, data=CAE) 

Les deux variables semblent être corrélées linéairement mais il faut le tester et le valider.

hist(CAE$tarse)

hist(CAE$poids)

Tester la normalite des variables.

shapiro.test(CAE$tarse)

Cette variable se distribue normalement (non rejet de H0 car la pvalue>0.05).

shapiro.test(CAE$poids) 

Cette variable, par contre, ne se distribue pas normalement (W = 0.96254, p-value = 0.002146). Comme l’une des deux variables ne suit pas une loi Normale, nous ne pouvons pas appliquer le test de correlation parametrique (corrélation linéaire) de Pearson.

Nous allons donc choisir de tester la corrélation non-paramétrique de Spearman ou de Kendall (non normalité des variables ou non linéarite de la covariation).

Le principe de calcul de ces corrélations non-paramétriques est le même que celui du calcul de Pearson, mais ici il est realisé sur les rangs, mais R fait tout seul le calcul des rangs

Tester la corrélation de ces variables avec cette écriture sans a priori sur le sens de la covariation.

cor.test(CAE$poids,CAE$tarse, method ="spearman")  

Utiliser l’argument alternative = “greater” pour indiquer un a priori sur le sens de la relation.

coeffSpearman<-cor.test(CAE$tarse,CAE$poids, method = "spearman",alternative ="greater") 
coeffSpearman

La covariation existe (HO rejeté car S = 135557, p-value = 2.759e-09), le sens est positif.

attributes(coeffSpearman)
coeffSpearman$estimate^2

Nous obtenons ainsi la valeur de l’intensité de la covariation (ici 0.2549623) entre le poids et le tarse des oisillons.