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

7.4.1 Présentation et donné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

 

7.4.2 Analyse des résultats

    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)

exp0741.jpg (55 Ko)

 

    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)
exp074$Bloc <- factor(exp074$Bloc)
exp074$Trait <- factor(exp074$Trait)
model <- aov(log10(Rend)~Var*Trait + Error(Bloc + Bloc:Var), data=exp074)
summary(model)


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"))
tapply(log10(exp074$Rend),exp074$Trait,mean)

       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)
cat("Efficacité", "\n",
"du test relatif aux variétés :", CMr/0.014073, "\n",
"du test relatif aux traitements et à l'interaction :", CMr/0.003101, "\n")

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)}
tapply(Rend,Trait,f) # coefficients de variation

        1         2         3 
0.2906801 0.3070779 0.3262702 
 
Haut de la page
Autres illustrations avec R





Dernière mise à jour : janvier 2005