Principes d'expérimentation :
planification des expériences et analyse de leurs résultats
Utilisation de R pour l'exemple du paragraphe 7.4
par Emmanuel Nowak
| 7.4.1 Présentation et données |
| 7.4.2 Analyse des résultats |
7.4 Exemple 1 : expérience en blocs aléatoires complets et parcelles divisées
1° Présentation
Se référer au paragraphe correspondant du livre. Il s'agit d'étudier les rendements, en kg de gousses fraîches par parcelle, de 6 variétés de haricots en fonction de 3 traitements appliqués à ces variétés. L'ensemble constitue une expérience factorielle complète comportant 18 objets.
2° Dispostif expérimental
On dispose de 72 petites parcelles, réparties en 4 blocs, constitués chacun de 6 grandes parcelles. C'est le principe d'une expérience en split-plot exposée au paragraphe 7.1.1 du livre.
3° Données
Il faut d'abord importer les données fournies dans le fichier 'exp074.txt' et diviser les valeurs des rendements qui y sont contenues par 100 pour obtenir les valeurs réelles en kg par parcelle. On suppose que ce fichier a été chargé dans 'C:\Dagnelie'. Par ailleurs, une colonne supplémentaire 'Objet' sera utile pour la suite :
| exp074 <-
read.table("C:/Dagnelie/expdonn/txt/exp074.txt",sep="\t",header=T) exp074$Rend <- exp074$Rend/100 exp074$Objet <- rep(1:18,each=4) exp074[1:5,] |
Var Trait Bloc Rend Objet 1 1 1 1 7.20 1 2 1 1 2 6.70 1 3 1 1 3 5.87 1 4 1 1 4 5.35 1 5 1 2 1 7.38 2 |
On peut calculer les rendements moyens en kg/parcelle pour chaque objet (4 répétitions) :
| Moy <- tapply(exp074$Rend,exp074$Objet,mean) cbind("Var"=rep(1:6,each=3),"Trait"=rep(1:3,length=18),"Moyennes"=Moy) |
Var Trait Moyennes 1 1 1 6.2800 2 1 2 7.3150 3 1 3 8.1600 4 2 1 6.8425 5 2 2 7.5250 6 2 3 8.0200 7 3 1 4.6250 8 3 2 6.2900 9 3 3 5.1500 10 4 1 4.8175 11 4 2 5.5350 12 4 3 6.0575 13 5 1 4.1900 14 5 2 4.3475 15 5 3 4.4275 16 6 1 3.3625 17 6 2 3.5875 18 6 3 3.5200 |
1° Examen préliminaire
Les représentations graphiques des données s'obtiennent à l'aide de la fonction 'stripchart', puisqu'on souhaite avoir le rendement (quantitatif) en fonction de la variété (qualitative) ou du traitement. La fonction 'plot' aurait été utilisable puisque les données sont codées numériquement, mais mal adaptée. La fonction 'par' permet de modifier certains paramètres graphiques, en particulier de diviser la fenêtre, ici en deux (l'agrandir ensuite).
| attach(exp074) par(mfrow=c(1,2),las=1) stripchart(Rend~Var, vertical=T, pch=21, xlab="Variétés", ylab="Rendements (kg/parc.)") stripchart(Rend~Trait, vertical=T, pch=21, xlab="Traitements", ylab="Rendements (kg/parc.)") detach(exp074) |
2° Analyse de la variance
Il faut ici réaliser une analyse de la variance à trois facteurs croisés, en définissant un terme correspondant à la partie aléatoire du modèle mixte du split-plot. Les interactions non spécifiées rentrent dans la composante résiduelle finale, appelée erreur 2 pour le split-plot. Notons qu'une transformation logarithmique doit être appliquée aux données initiales afin de stabiliser les variances.
exp074$Var <- factor(exp074$Var) |
Error: Bloc
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 3 0.147833 0.049278
Error: Bloc:Var
Df Sum Sq Mean Sq F value Pr(>F)
Var 5 1.05529 0.21106 14.997 2.238e-05 ***
Residuals 15 0.21110 0.01407
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Trait 2 0.056757 0.028378 9.1505 0.0006122 ***
Var:Trait 10 0.044062 0.004406 1.4208 0.2107896
Residuals 36 0.111646 0.003101
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
|
Le carré moyen du facteur blocs vaut 0,0493 et ne sert de base de comparaison pour aucun test. Le carré moyen de l'interaction blocs-variétés vaut 0,0141 et sert de base de comparaison pour tester l'effet du facteur variétés. Le carré moyen des termes résiduels (within error) vaut 0,0031 et sert de base de comparaison pour tester les effets du facteur traitements et de l'interaction variétés-traitements.
Remarque : avant d'interpréter les résultats de l'analyse de la variance, il faut en vérifier la validité des conditions d'application. La procédure est exactement la même que celle de l'exemple du paragraphe 5.4. Elle n'est donc pas reprise dans cette page.
3° Différence de moyennes
Les moyennes des logarithmes des 24 observations relatives à chacun des traitements valent :
exp074$Trait <- factor(exp074$Trait,
labels=c("témoin","fumure azotée","inoculation")) |
témoin fumure azotée inoculation
0.6814398 0.7387542 0.7430149
|
Afin de déterminer les intervalles de confiance pour les différences entre témoin et inoculation puis témoin et fumure azotée, il ne manque que le quantile de la loi de Student avec 36 degrés de liberté :
qt(0.025,df=36,lower.tail=F) |
[1] 2.028094 |
Les intervalles de confiance pour l' influence de l'inoculation et celle de la fumure azotée sont donc centrés respectivement sur 0,73875-0,68144 = 0,05731 et 0,74301-0,68144 = 0,06157 et leur demi-amplitude vaut :
2.028*sqrt(2*0.003101/24) |
[1] 0.03260080 |
Les intervalles de confiances respectifs sont alors [0,02471 ; 0,08991] et [0,02897 ; 0,09417]. Rappelons que l'on travaille ici sur les données transformées. On peut en déduire des intervalles de confiance pour les rapports des rendements moyens en utilisant la fonction inverse du logarithme décimal :
10^c(0.02471,0.08991) # inoculation / témoin |
[1] 1.058547 1.230014 |
10^c(0.02897,0.09417) # fumure azotée / témoin |
[1] 1.068981 1.242138 |
Notons que cette dernière opération, dont la légitimité est loin d'être évidente, nécessite quelques explications qui seront données au paragraphe 5°.
4° Efficacité relative
Les formules du paragraphe 6.3 mènent à :
CMr <-
(0.212202+0.111646)/(15+36) |
Efficacité du test relatif aux variétés : 0.4512159 du test relatif aux traitements et à l'interaction : 2.047714 |
5° À propos de la transformation logarithmique
Les rendements moyens que l'on peut obtenir après transformation logarithmique puis retour à la variable initiale sont en fait des moyennes géométriques et non arithmétiques. En ce qui concerne les intervalles de confiance des différences, en notant X1 et X2 les variables rendements avec chacun des deux traitements considérés et E(.) l'espérance mathématique, on obtient tout d'abord un intervalle de confiance pour
E(log(X2)) - E(log(X1)).
Or on cherche un intervalle de confiance pour
log(E(X2)) - log(E(X1)),
ce qui n'est pas la même chose. Cependant, un développement en série de Taylor de log(X) autour de E(X), limité aux trois premiers termes montre que approximativement :
E(log(X)) = log(E(X)) - CV2/2.loge(10)
où CV est le coefficient de variation de X et loge le logarithme népérien. Si les coefficients de variations sont semblables avec les deux traitements, alors les biais se compensent et on peut approximer log(E(X2))-log(E(X1)) par E(log(X2))-E(log(X1)).
f <- function(x){sd(x)/mean(x)} |
1 2 3 0.2906801 0.3070779 0.3262702 |
| Haut de la page |
| Autres illustrations avec R |
Dernière mise à jour : janvier 2005