Principes d'expérimentation :

planification des expériences et analyse de leurs résultats

 
Utilisation de R pour l'exemple du paragraphe 8.6

par Emmanuel Nowak
8.6.1 Présentation et données
8.6.2 Analyse des résultats

 

 

8.6 Exemple 2 : expérience avec permutation des objets (cross-over)

8.6.1 Présentation et données

    1° Présentation

Se référer au paragraphe correspondant du livre. Il s'agit de comparer trois alimentations A, B et C dont la première est une alimentation témoin. L'expérience a été réalisée selon un dispositif en cross-over faisant intervenir 6 vaches et 3 périodes d'observation, chaque alimentation étant administrée une fois à chaque animal et deux fois pour chaque période. Le dispositif est équilibré pour les effets résiduels.

    2° Données

Il faut d'abord importer les données fournies dans le fichier 'exp086.txt' et diviser les valeurs des productions qui y sont contenues par 10 pour obtenir les valeurs réelles en kg/jour. On suppose que ce fichier a été chargé dans 'C:\Dagnelie'. Par ailleurs, la variable 'Alim' sera recodée afin de conserver les notations du livre.

exp086 <- read.table("C:/Dagnelie/expdonn/txt/exp086.txt",sep="\t",header=T)
exp086$Prod <- exp086$Prod/10
exp086$Alim <- factor(exp086$Alim, labels=c("A","B","C"))
exp086[1:5,]

  Alim Per Vache Prod Resid
1    B   1     1 21.7     0
2    C   1     2 20.4     0
3    A   1     3 25.9     0
4    B   1     4 23.5     0
5    A   1     5 19.1     0

 

8.6.2 Analyse des résultats

    1° Examen préliminaire

La représentation graphique de la production en fonction de l'alimentation s'obtient à l'aide de la fonction 'stripchart'.

attach(exp086)
par(las=1)
stripchart(Prod~Alim, ylim=c(15,30), vert=T, pch=21, xlab="Alimentations", ylab="Productions (kg/jour)")

exp0861.jpg (21,2 Ko)

 

La procédure est la même pour la production en fonction des deux autres facteurs. La fonction 'par' permet de modifier certains paramètres graphiques, en particulier de diviser la fenêtre, ici en deux (l'agrandir ensuite).

par(mfrow=c(1,2),las=1)
stripchart(Prod~Per, ylim=c(15,30), vert=T, pch=21, xlab="Périodes", ylab="Productions (kg/jour)")
stripchart(Prod~Vache, ylim=c(15,30), vert=T, pch=21, xlab="Vaches", ylab="Productions (kg/jour)")

exp0862.jpg (42,5 Ko)

 

    2° Analyse de la variance et test d'additivité

Il faut ici réaliser une analyse de la variance à trois facteurs croisés, en supposant que le modèle est additif puisque l'effet du facteur périodes ne peut être considéré comme étant aléatoire. Les variables périodes et vaches étant codées numériquement, il ne faut pas oublier de préciser de les considérer comme des facteurs.

model1 <- aov(Prod~Alim+factor(Per)+factor(Vache))
summary(model1)

              Df Sum Sq Mean Sq F value    Pr(>F)    
Alim           2 12.988   6.494  10.014 0.0066370 ** 
factor(Per)    2 13.778   6.889  10.623 0.0055984 ** 
factor(Vache)  5 71.698  14.340  22.113 0.0001721 ***
Residuals      8  5.188   0.648                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

L'analyse de la variance indique une influence très significative de chacun des facteurs pris en compte.

La condition d'additivité du modèle peut être contrôlée à l'aide du test de Tukey. Il faut pour cela calculer la composante de non-additivité  puis la valeur du F de Fisher qui s'en déduit, avec 1 et  8-1 = 7 degrés de liberté :

SCEres <- 17*var(model1$res) # qui vaut 5,188 en arrondissant
y <- (Prod-model1$res)^2
model2 <- aov(y~Alim+factor(Per)+factor(Vache))
SCEadd <- sum(y*model1$res)^2 / (17*var(model2$res))
F <- SCEadd/((SCEres-SCEadd)/7)
p <- pf(F, df1=1, df2=7, lower.tail=F)
cat("F=",F," p=",p, "\n")

F= 0.4320448  p= 0.5319969

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° Efficacité relative

En considérant l'expérience comme étant complètement aléatoire, les sommes des carrés relatives aux facteurs périodes et vaches s'ajoutent à la résiduelle, d'où l'efficacité relative :

((13.7778+71.6978+5.1878)/(2+5+8))/0.6485

[1] 9.320319

 

    4° Interprétation

Le facteur alimentation étant hautement significatif, la question est maintenant de comparer les alimentations B et C à l'alimentation A. S'agissant de comparaisons unilatérales de plusieurs moyennes à un témoin, dans l'optique d'identifier la ou les productions supérieures à celle du témoin, la méthode à utiliser est celle de Dunnett. Les moyennes observées sont :

tapply(Prod,Alim,mean)

       A        B        C 
20.58333 21.30000 19.25000

La valeur de t intervenant dans le calcul de la plus petite différence significative au sens de Dunnett, lue dans la table correspondant à des comparaisons unilatérales, vaut 2,22 pour trois populations et huit degrés de liberté. Cette plus petite différence est alors égale à :

round(2.22*sqrt(2*0.648/6), digit=2)

[1] 1.03 

Aucune des deux alimentations B et C n'est donc significativement supérieure à l'alimentation A.

 

    5° Effets résiduels

L'ajout du facteur effets résiduels permet de tester l'existence éventuelle de tels effets.

model3 <- aov(Prod~Alim+factor(Per)+factor(Vache)+factor(Resid))
summary(model3)

              Df Sum Sq Mean Sq F value    Pr(>F)    
Alim           2 12.988   6.494 16.6926 0.0035355 ** 
factor(Per)    2 13.778   6.889 17.7080 0.0030405 ** 
factor(Vache)  5 71.698  14.340 36.8600 0.0001984 ***
factor(Resid)  2  2.854   1.427  3.6676 0.0910860 .  
Residuals      6  2.334   0.389                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

Les effets résiduels n'apparaissent donc pas comme significatifs.

 

    6° Modèles des analyses de la variance

Notons ici l'importance de l'ordre d'apparition des facteurs dans le dernier modèle d'analyse de la variance, celle-ci étant non orthogonale. Le facteur effets résiduels doit y intervenir en dernier lieu, sous peine d'obtenir une fausse décomposition de la somme des carrés des écarts, comme l'illustre l'exemple ci-dessous.

model4 <- aov(Prod~factor(Resid)+Alim+factor(Per)+factor(Vache)) # mauvais ordre !
summary(model4)

              Df Sum Sq Mean Sq F value    Pr(>F)    
factor(Resid)  3 15.486   5.162 13.2691 0.0046677 ** 
Alim           2 13.347   6.674 17.1543 0.0032981 ** 
factor(Per)    1  0.333   0.333  0.8568 0.3903378    
factor(Vache)  5 72.150  14.430 37.0927 0.0001949 ***
Residuals      6  2.334   0.389                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
 
Haut de la page
Autres illustrations avec R





Dernière mise à jour : janvier 2005