Principes d'expérimentation :

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

 
Utilisation de R pour l'exemple du paragraphe 5.4

par Emmanuel Nowak
5.4.1 Présentation et données
5.4.2 Analyse des résultats
5.4.3 Importance de la randomisation

 

 

5.4 Exemple 1 : expérience complètement aléatoire à deux facteurs

5.4.1 Présentation et données

    1° Présentation

Se référer au paragraphe correspondant du livre. Il s'agit d'étudier l'influence de la dimension et du degré d'humidité de fragments de bois soumis à la carbonisation sur le rendement en charbon de bois.

    2° Données

Il faut d'abord importer les données fournies dans le fichier 'exp054.txt'. On suppose que ce fichier a été chargé dans 'C:\Dagnelie'. Il reste ensuite à diviser les valeurs des rendements qui y sont contenues par 100 pour obtenir les valeurs réelles :

exp054 <- read.table("C:/Dagnelie/expdonn/txt/exp054.txt",sep="\t",header=T)
exp054$Rend <- exp054$Rend/100
exp054[1:5,]

  Dimen Humid Repet  Rend
1     2     0     1 30.00
2     2     0     2 29.67
3     2     0     3 29.78
4     2    10     1 29.82
5     2    10     2 29.71 

 

5.4.2 Analyse des résultats

    1° Examen préliminaire

La représentation graphique des données s'obtient en utilisant la fonction 'plot'. La fonction 'par' permet de modifier certains paramètres graphiques, en particulier de diviser la fenêtre graphique en plusieurs parties, ici en deux (l'agrandir ensuite).

attach(exp054)
par(mfrow=c(1,2))
plot(Dimen,Rend,xlab="Dimensions (cm)", ylab="Rendements (%)",las=1)
plot(Humid,Rend,xlab="Humidités (%)", ylab="Rendements (%)",las=1)

exp0541.jpg (47 Ko)

La valeur 33,11 doit être testée par la méthode de Grubbs. Il faut d’abord effectuer l’analyse de variance à deux critères à l’aide de la fonction 'aov', afin d’en calculer les résidus réduits. Les variables 'Dimen' et Humid' étant continues, il faut préciser de les considérer comme des facteurs qualitatifs :

model1 <- aov(Rend~as.factor(Dimen)*as.factor(Humid))
summary(model1)


                                  Df  Sum Sq Mean Sq F value  Pr(>F)  
as.factor(Dimen)                   2  3.6184  1.8092  3.2509 0.05632 .
as.factor(Humid)                   3  1.0922  0.3641  0.6542 0.58818  
as.factor(Dimen):as.factor(Humid)  6  1.4771  0.2462  0.4424 0.84309  
Residuals                         24 13.3566  0.5565                  
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

La variance résiduelle est estimée à 0,5565. Il ne reste plus qu’à calculer le 10ème résidu réduit, correspondant à la valeur suspecte :

model1$res[10]/sqrt(0.5565) # à comparer à 2.878 lue dans la table
detach(exp054)


      10
3.074218

 
La valeur est donc considérée comme anormale et ne sera plus prise en considération. Le fichier sans cette observation sera nommé 'exp054b' :

exp054b <- exp054[-10,]
attach(exp054b)

par(mfrow=c(1,2))
plot(Dimen,Rend,xlab="Dimensions (cm)", ylab="Rendements (%)",las=1)
plot(Humid,Rend,xlab="Humidités (%)", ylab="Rendements (%)",las=1)

detach(exp054b)

exp0542.jpg (53 Ko)

De nombreuses options existent pour les graphiques : consulter éventuellement l'aide sur 'plot' ou 'plot.default' ou encore 'par'.


    2° Analyse de la variance

La valeur considérée comme manquante doit être estimée par la moyenne des autres observations relatives au même objet, soit 29,67. Le calcul est immédiat, mais si les répétitions avaient été plus nombreuses, on aurait procédé de la façon suivante :

exp054c <- exp054
exp054c$Rend[10] <- mean(exp054$Rend[exp054$Dimen==2 & exp054$Humid==40 & exp054$Rep!=1])

Le fichier corrigé est appelé 'exp054c' et il ne reste plus qu’à effectuer l’analyse de la variance :

attach(exp054c)
model2 <- aov(Rend~as.factor(Dimen)*as.factor(Humid))
summary(model2)
                                  Df Sum Sq Mean Sq F value  Pr(>F)  
as.factor(Dimen)                   2 1.2925  0.6463  2.8369 0.07836 .
as.factor(Humid)                   3 0.4367  0.1456  0.6389 0.59735  
as.factor(Dimen):as.factor(Humid)  6 0.1947  0.0325  0.1425 0.98886  
Residuals                         24 5.4675  0.2278                  
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Attention : les sommes des carrés sont justes mais il faut recalculer le carré moyen résiduel qui n’a que 23 degrés de liberté et non 24 du fait de l’estimation de la valeur manquante et refaire les calculs :

model2$df.residual <- 23
summary(model2)
                                  Df Sum Sq Mean Sq F value  Pr(>F)  
as.factor(Dimen)                   2 1.2925  0.6463  2.7187 0.08713 .
as.factor(Humid)                   3 0.4367  0.1456  0.6123 0.61392  
as.factor(Dimen):as.factor(Humid)  6 0.1947  0.0325  0.1365 0.98998  
Residuals                         23 5.4675  0.2377                  
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
    

 

    3° Conditions d'application

Afin de vérifier la normalité de la distribution des résidus, une représentation de ceux-ci sous forme de box-plot peut être obtenue à l'aide de la fonction 'boxplot'.

 boxplot(model2$res,horizontal=T,xlab="Résidus (%)")

exp0543.jpg (18 Ko)

Le diagramme de probabilité s'obtient quant à lui à l'aide de la fonction 'qqnorm', cette dernière représentant les quantiles observés relatifs aux résidus en fonction des quantiles théoriques aussi appelés scores normaux, et non l'inverse. Pour obtenir les scores normaux en fonction des résidus, il suffit de permuter les axes.

par(mfrow=c(2,2))
boxplot(model2$res,horizontal=T,xlab="Résidus(%)")
diagramme <- qqnorm(model2$res,main="") 
# permutation des axes :
plot(diagramme$y,diagramme$x, xlab="Résidus (%)", ylab="Scores normaux")

exp0544.jpg (18 Ko)


Afin de vérifier l'homogénéité des variances, le test de Levene consiste en une analyse de la variance sur les valeurs absolues des résidus :

model3 <- aov(abs(model2$res)~as.factor(Dimen)*as.factor(Humid))
summary(model3)
                                  Df  Sum Sq Mean Sq F value Pr(>F)
as.factor(Dimen)                   2 0.08225 0.04113  0.8158 0.4542
as.factor(Humid)                   3 0.17556 0.05852  1.1609 0.3452
as.factor(Dimen):as.factor(Humid)  6 0.48632 0.08105  1.6079 0.1882
Residuals                         24 1.20986 0.05041          

Il reste à refaire les calculs, comme précédemment, en utilisant 23 degrés de liberté pour le carré moyen résiduel.

model3$df.residual <- 23
summary(model3)
                                  Df  Sum Sq Mean Sq F value Pr(>F)
as.factor(Dimen)                   2 0.08225 0.04113  0.7818 0.4693
as.factor(Humid)                   3 0.17556 0.05852  1.1125 0.3644
as.factor(Dimen):as.factor(Humid)  6 0.48632 0.08105  1.5409 0.2096
Residuals                         23 1.20986 0.05260 

 

    4° Interprétation

L'analyse de la variance ne fait apparaître aucune différence significative, mais l'approche par régression conduit à des conclusions différentes en ce qui concerne l'influence de la dimension des fragments de bois :

model4 <- lm(Rend~Dimen)
anova(model4)
Analysis of Variance Table

Response: Rend
          Df Sum Sq Mean Sq F value  Pr(>F)  
Dimen      1 1.1978  1.1978  6.5752 0.01493 *
Residuals 34 6.1937  0.1822                  
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 


La somme des carrés relative à la relation linéaire vaut 1,1978. Elle est égale au carré moyen correspondant (un seul degré de liberté), à diviser par le carré moyen résiduel (=0,2377) de l'analyse de la variance pour obtenir F=5,04. Le degré de signification vaut alors :

pf(5.04,df1=1,df2=23,lower.tail=F)
[1] 0.03468056

Les coefficients de l’équation de régression linéaire sont  contenus dans l’objet 'model4' :

model4$coef
detach(exp054c)
(Intercept)       Dimen 
  29.917083   -0.073125 

 

Haut de la page

5.4.3 Importance de la randomisation

    1° Principe

Les données qui auraient été obtenues en présence des dérives linéaires envisagées doivent être calculées en partant du fichier initial 'exp054'. On ajoute tout d'abord les dérives, positives ou négatives, aux valeurs des rendements de départ :

exp054d <- exp054
exp054d$Neg.01 <- exp054$Rend-seq(0,35*0.01, by=0.01)
exp054d$Neg.005 <- exp054$Rend-seq(0,35*0.005, by=0.005)
exp054d$Pos.005 <- exp054$Rend+seq(0,35*0.005, by=0.005)
exp054d$Pos.01 <- exp054$Rend+seq(0,35*0.01, by=0.01)
exp054d[1:6,]
  Dimen Humid Repet  Rend Neg.01 Neg.005 Pos.005 Pos.01
1     2     0     1 30.00  30.00  30.000  30.000  30.00
2     2     0     2 29.67  29.66  29.665  29.675  29.68
3     2     0     3 29.78  29.76  29.770  29.790  29.80
4     2    10     1 29.82  29.79  29.805  29.835  29.85
5     2    10     2 29.71  29.67  29.690  29.730  29.75
6     2    10     3 29.87  29.82  29.845  29.895  29.92

Puis il convient d'éliminer la 10ème observation, considérée comme aberrante, et de l'estimer par la moyenne des autres observations relatives au même objet (dimension = 2 et humidité = 40) :

exp054d[10,4:8] <- mean(exp054d[exp054d$Dimen==2 & exp054d$Humid==40 & exp054d$Rep!=1,4:8])

 

    2° Résultats

Il faut maintenant effectuer les quatre analyses de la variance relatives aux données modifiées par les différentes dérives linéaires afin d'obtenir les sommes des carrés des écarts, puis enlever un degré de liberté à la composante résiduelle car une valeur a été estimée et recalculer les F-values ainsi que les p-values. On peut aussi calculer les degrés de signification relatifs à la composante linéaire du facteur dimension.

Les résultats obtenus ici peuvent différer très légèrement de ceux du livre, dans lequel la valeur du rendement estimé a été arrondie à deux ou trois décimales, soit 29,56,  29,618,  29,722 et 29,78. La procédure étant purement répétitive et la même que celle utilisée précédemment, l'ensemble des calculs n'est pas repris dans cette page.

Haut de la page
Autres illustrations avec R

Voir les calculs



Dernière mise à jour : janvier 2005