Statistique théorique et appliquée - Tome 2
Utilisation de R pour les exemples 9.2.1 à 9.3.2
| Exemples 9.2.1 à 9.3.1 |
| Exemple 9.3.2 |
par Emmanuel Nowak
Il s'agit ici de comparer les hauteurs moyennes des arbres de trois types de hêtraies. Il faut tout d'abord importer les données relatives à l'exemple 2.3.1, fournies dans le fichier 's2e02031.txt', diviser les hauteurs par 10 afin d'obtenir des valeurs en mètres et déclarer la variable 'Types' comme facteur :
| s2e02031 <- read.table("C:/Dagnelie/st2donn/txt.2/s2e02031.txt",sep="\t",header=T) s2e02031$Haut <- s2e02031$Haut/10 s2e02031$Types <- factor(s2e02031$Types) summary(s2e02031) |
Types Haut
1:13 Min. :18.90
2:14 1st Qu.:23.70
3:10 Median :25.00
Mean :24.98
3rd Qu.:26.70
Max. :28.50
|
La moyenne générale, indiquée ci-dessus parmi les statistiques élémentaires, est de 24,98 m. On peut facilement obtenir les moyennes par types de hêtraies :
| attach(s2e02031) tapply(Haut,Types,mean) |
1 2 3 25.96923 25.38571 23.14000 |
Le tableau d'analyse de la variance, contenant les éléments de décomposition de la variation totale ainsi que les aspects inférentiels, permet de comparer les hauteurs moyennes. Il s'agit ici d'un modèle fixe d'analyse de la variance à un critère de classification :
| model <- aov(Haut~Types) summary(model) |
Df Sum Sq Mean Sq F value Pr(>F) Types 2 48.881 24.441 7.1238 0.002606 ** Residuals 34 116.649 3.431 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 |
Le logiciel R permet également d'obtenir les limites de confiance des différences en utilisant la méthode HSD de Tukey. Cette méthode, plus élaborée que celle de Bonferroni évoquée dans le livre, tient compte du problème de comparaisons multiples. Les intervalles de confiance, que l'on a choisi d'ordonner ci-dessous, sont alors plus étendus que ceux obtenus sans correction :
| TukeyHSD(model,ordered=T) |
Tukey multiple comparisons of means
95% family-wise confidence level
factor levels have been ordered
Fit: aov(formula = Haut ~ Types)
$Types
diff lwr upr
2-3 2.2457143 0.3664596 4.124969
1-3 2.8292308 0.9200973 4.738364
1-2 0.5835165 -1.1646799 2.331713
|
L'utilisation de l'analyse de la variance impose de vérifier que les conditions de normalité et d'égalité des variances sont respectées, et d'identifier d'éventuelles observations aberrantes. Plusieurs méthodes peuvent être mises en oeuvre dans ce but et ont fait l'objet des exemples 3.3.1, 3.3.3, 3.5.2 et 7.5.1. Les données sont toujours celles de l'exemple 2.3.1 et ont été importées précédemment.
La normalité peut être contrôlée par la méthode des diagrammes de probabilité (exemple 3.3.1). Il faut tout d'abord définir un tableau contenant les hauteurs et les quantiles normaux associés, ceci pour chaque type de hêtraie :
| diag.prob <- tapply(Haut,Types,qqnorm) names(diag.prob) <- c("type.1","type.2","type.3") diag.prob$type.3 # valeurs pour le type 3 |
$x [1] -1.5466353 -1.0004905 -0.6554235 -0.3754618 -0.1225808 0.1225808 0.3754618 0.6554235 1.0004905 1.5466353 $y [1] 18.9 21.1 21.2 22.1 22.5 23.6 24.5 24.6 26.2 26.7 |
On voit ci-dessus que les hauteurs sont en y et les quantiles en x. Les axes des diagrammes de probabilités seront donc inversés par rapport à ceux du livre :
| plot(diag.prob$type.1,xlim=c(-2,2),ylim=c(18,29),pch=21,bg="red",xlab="Quantiles normaux",ylab="Hauteurs (m)") points(diag.prob$type.2,pch=21,bg="green") points(diag.prob$type.3,pch=21,bg="blue") legend(1,21,c("Type 1","Type 2","Type 3"),pt.bg=c("red","green","blue"), pch=21) |
On peut également vérifier la normalité, pour chaque type de hêtraie, à l'aide du test de Shapiro et Wilk (exemple 3.3.3). On obtient par exemple, pour le type 3 :
| shapiro.test(Haut[Types=="3"]) |
Shapiro-Wilk normality test data: Haut[Types == "3"] W = 0.9722, p-value = 0.9104 |
Il est également conseillé d'identifier d'éventuelles observations aberrantes à l'aide du test de Grubbs (exemple 3.5.2). Le carré moyen résiduel de l'analyse de la variance peut servir d'estimation de la variance, afin de déterminer la valeur maximale pour les valeurs absolues des résidus réduits :
| max(abs(model$res)/sqrt(3.431)) |
[1] 2.289051 |
Afin de retrouver les résidus réduits du livre, il suffit de remplacer le carré moyen résiduel par la variance des résidus. Ceci est à comparer à la valeur limite :
| qt(1-0.05/74,df=34)*sqrt(34/(33+qt(1-0.05/74,df=34)^2)) |
3.028634 |
Enfin, il reste à vérifier l'égalité des variances. Le test de Levene consiste à réaliser une analyse de la variance sur les valeurs absolues des résidus :
| summary(aov(abs(model$res)~Types)) |
Df Sum Sq Mean Sq F value Pr(>F) Types 2 3.8003 1.9002 2.2262 0.1234 Residuals 34 29.0206 0.8535 |
Un test de Bartlett aurait conduit à la même conclusion concernant les variances :
| bartlett.test(model$res~Types) |
Bartlett test for homogeneity of variances data: model$res by Types Bartlett's K-squared = 3.4657, df = 2, p-value = 0.1768 |
Avant de traiter d'autres exemples, ne pas oublier de détacher le fichier de travail :
| detach(s2e08051) |
| Haut de la page |
| Autres illustrations avec R |
Dernière mise à jour : août 2006