Principes d'expérimentation :

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

 
Utilisation de R pour l'exemple du paragraphe 6.5

par Emmanuel Nowak
6.5.1 Présentation et données
6.5.2 Analyse des résultats : analyse de la variance
6.5.3 Analyse des résultats : régression

 

 

6.5 Exemple 2 : expérience en blocs aléatoires complets

6.5.1 Présentation et données

    1° Présentation

Se référer au paragraphe correspondant du livre. Il s'agit de comparer des fumures appliquées sur blé au Rwanda, dans des conditions de fertilité naturelle extrêmement faible. L'analyse statistique porte sur le rendement en fonction des niveaux de phosphore (100, 200 et 300 kg/ha) et de calcium (1000, 4500 et 8000 kg/ha). Il y a donc 9 objets, répétés sur 3 blocs.

    2° Données

Il faut d'abord importer les données fournies dans le fichier 'exp065.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' est à définir pour la numérotation des objets, de 8 à 16 pour reprendre les notations du livre :

exp065 <- read.table("C:/Dagnelie/expdonn/txt/exp065.txt",sep="\t",header=T,na.strings="*")
exp065$Rend <- exp065$Rend/100
exp065$Objet <- rep(8:16,each=3)
exp065[1:5,]

  Phos Calc Bloc Rend Objet
1  100 1000    1 1.30     8
2  100 1000    2 0.80     8
3  100 1000    3 2.01     8
4  100 4500    1 2.05     9
5  100 4500    2 2.37     9

On peut calculer les rendements moyens en kg/parcelle pour chaque objet (3 répétitions) ou en t/ha (il suffit de diviser le rendement en kg/parcelle par 1,8) :

Moy <- tapply(exp065$Rend, exp065$Objet, mean,na.rm=T)
Moy <- cbind("Moyennes (kg/p)"=Moy, "Moyennes (t/ha)"=Moy/1.8)
rownames(Moy) <- paste("Objet",8:16)
round(Moy,digit=3)

         Moyennes(kg/p) Moyennes(t/ha)
Objet 8           1.370          0.761
Objet 9           2.313          1.285
Objet 10          2.307          1.281
Objet 11          2.535          1.408
Objet 12          3.180          1.767
Objet 13          2.980          1.656
Objet 14          2.653          1.474
Objet 15          3.333          1.852
Objet 16          3.390          1.883

 

6.5.2 Analyse des résultats : analyse de la variance

    1° Examen préliminaire

La représentation graphique des données s'obtient à l'aide de la fonction 'plot', utilisée ici avec quelques options. Le premier graphique représente le rendement en fonction de la quantité de phosphore, pour chacun des trois niveaux de calcium, toutes les données étant converties en t/ha :

attach(exp065)

plot(Phos/1000,Rend/1.8, axes=F, frame.plot=T,
xlab="Phosphore (t/ha)",ylab="Rendements(t/ha)",
ylim=c(0.4,2.2),pch=21,bg=c("red", "green", "blue")[as.numeric(as.factor(Calc))])

axis(1, at=c(0.1,0.2,0.3))
axis(2, at=seq(0.5,2,by=0.5), las=1)
legend(0.22,0.92,c("Ca=1.0 t/ha","Ca=4.5 t/ha","Ca=8.0 t/ha"), pt.bg=c("red","green","blue"), pch=21)

exp0651.jpg (35 Ko)

 

Le second graphique représente le rendement en fonction de la quantité de calcium, pour chacun des trois niveaux de phosphore, toutes les données étant converties en t/ha :

plot(Calc/1000,Rend/1.8, axes=F, frame.plot=T,
xlab="Calcium (t/ha)",ylab="Rendements(t/ha)",
ylim=c(0.4,2.2),pch=21,bg=c("red", "green", "blue")[as.numeric(as.factor(Phos))])

axis(1, at=c(1,4.5,8))
axis(2, at=seq(0.5,2,by=0.5), las=1)
legend(4.72,0.92,c("P2O5=0.1 t/ha","P2O5=0.2 t/ha","P2O5=0.3 t/ha"),
pt.bg=c("red","green","blue"), pch=21)

detach(exp065)
 

exp0652.jpg (38 Ko)



    2° Analyse de la variance

Il faut d'abord estimer la donnée manquante relative à l'objet 11, bloc 1. On peut faire facilement le calcul à la main mais la procédure générale est la suivante :

p1 <- length(table(exp065$Objet))
p2 <- length(table(exp065$Bloc))
somme1 <- sum(na.omit(exp065$Rend[exp065$Objet==11]))
somme2 <- sum(na.omit(exp065$Rend[exp065$Bloc==1]))
somme <- sum(na.omit(exp065$Rend))

exp065$Rend[10] <- round( (p1*somme1 + p2*somme2 - somme) / ((p1-1)*(p2-1)), digit=2)
exp065$Rend[10]
 

[1] 2.25

 

La fonction 'aov' permet de réaliser l'analyse de la variance correspondant au modèle croisé mixte à trois critères de classification, et la fonction 'summary' présente les résultats sous la forme d'un tableau d'analyse de la variance :

attach(exp065)
model1 <- aov(Rend~as.factor(Phos)*as.factor(Calc)+as.factor(Bloc))
summary(model1)
 
 

                                Df Sum Sq Mean Sq F value    Pr(>F)    
as.factor(Phos)                  2 6.2949  3.1475 42.6080 3.899e-07 ***
as.factor(Calc)                  2 3.5022  1.7511 23.7052 1.643e-05 ***
as.factor(Bloc)                  2 1.9191  0.9596 12.9899 0.0004453 ***
as.factor(Phos):as.factor(Calc)  4 0.1525  0.0381  0.5162 0.7249848    
Residuals                       16 1.1819  0.0739                      
---
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 15 degrés de liberté et non 16, du fait de l’estimation de la valeur manquante, ainsi que les valeurs des variables de Fisher-Snedecor (F-values) et les degrés de signification associés Pr(>F) ou  p-values :

model1$df.residual <- 15
summary(model1)
 
 

                                 Df Sum Sq Mean Sq F value    Pr(>F)    
as.factor(Phos)                  2 6.2949  3.1475 39.9450 9.807e-07 ***
as.factor(Calc)                  2 3.5022  1.7511 22.2237 3.271e-05 ***
as.factor(Bloc)                  2 1.9191  0.9596 12.1781 0.0007213 ***
as.factor(Phos):as.factor(Calc)  4 0.1525  0.0381  0.4839 0.7473482    
Residuals                       15 1.1819  0.0788                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Remarques : 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. Par ailleurs, le facteur bloc ne devrait pas être testé par rapport à la variation dite résiduelle du modèle, qui regroupe en réalité trois interactions.

 

    3° Interprétation

Les deux facteurs pris en considération (phosphore et calcium) sont significatifs (degrés de signification proches de zéro) et leur interaction est négligeable.

Par ailleurs, l'efficacité relative du dispositif en blocs aléatoires complets par rapport à une répartition aléatoire, correspondant à une analyse de la variance à deux facteurs avec interaction, s'obtient en première approximation par un simple rapport des carrés moyens résiduels :

round(((1.9191+1.1819)/(2+15))/CMR ,digit=3)
[1] 2.315

Cependant, cette valeur est surestimée. Un calcul exact basé sur les relations du paragraphe 6.3.3° conduirait à une efficacité proche de 190%.

Voir les calculs

 

    4° Analyse complémentaire

La décomposition des sommes des carrés faisant intervenir les parties linéaires et quadratiques des facteurs 'Phosphore' et 'Calcium'  peut être obtenue par deux régressions simples ou par régression multiple. La fonction 'lm' permet d'ajuster le modèle linéaire et la fonction 'anova' présente les résultats sous la forme d'un tableau d'analyse de la variance :

model2 <- lm(Rend~Phos+I(Phos^2)+Calc+I(Calc^2)+as.factor(Phos):as.factor(Calc)+as.factor(Bloc))
anova(model2)
 
 

Analysis of Variance Table

Response: Rend
                                Df Sum Sq Mean Sq F value    Pr(>F)    
Phos                             1 5.7348  5.7348 77.6327 1.551e-07 ***
I(Phos^2)                        1 0.5602  0.5602  7.5834 0.0141236 *  
Calc                             1 2.4494  2.4494 33.1584 2.936e-05 ***
I(Calc^2)                        1 1.0528  1.0528 14.2521 0.0016572 ** 
as.factor(Bloc)                  2 1.9191  0.9596 12.9899 0.0004453 ***
as.factor(Phos):as.factor(Calc)  4 0.1525  0.0381  0.5162 0.7249848    
Residuals                       16 1.1819  0.0739                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

 

Attention : il faut refaire les calculs pour les valeurs de F et de p puisque le carré moyen résiduel n'a que 15 degrés de liberté, à cause de l'estimation de la donnée manquante.

model2$df.residual <- 15
anova(model2)
 
 

Analysis of Variance Table

Response: Rend
                                Df Sum Sq Mean Sq F value    Pr(>F)    
Phos                             1 5.7348  5.7348 72.7806 3.867e-07 ***
I(Phos^2)                        1 0.5602  0.5602  7.1094 0.0176065 *  
Calc                             1 2.4494  2.4494 31.0860 5.302e-05 ***
I(Calc^2)                        1 1.0528  1.0528 13.3613 0.0023442 ** 
as.factor(Bloc)                  2 1.9191  0.9596 12.1781 0.0007213 ***
as.factor(Phos):as.factor(Calc)  4 0.1525  0.0381  0.4839 0.7473482    
Residuals                       15 1.1819  0.0788                      
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Haut de la page

6.5.3 Analyse des résultats : régression

    1° Deux facteurs considérés séparément

Deux régressions multiples successives vont permettre d'ajuster des équations du deuxième degré, l'une pour le phosphore et l'autre pour le calcium. Les coefficients de régression seront contenus dans le vecteur 'coefficient' attaché au modèle correspondant. Afin de travailler avec des variables exprimées en t/ha, il faut diviser le rendement par 1,8 ainsi que les quantités de phosphore et de calcium par 1000. Ceci explique les légères différences entre les résultats obtenus ici et ceux du livre, dans lequel ces valeurs en t/ha ont été arrondies.

Pour le phosphore :

model3 <- lm(I(Rend/1.8)~I(Phos/1000)+I((Phos/1000)^2))
model3$coefficient
     (Intercept)     I(Phos/1000) I((Phos/1000)^2) 
       0.2864198        9.9259259      -16.9753086 

Pour le calcium :

model4 <- lm(I(Rend/1.8)~I(Calc/1000)+I((Calc/1000)^2))
model4$coefficient
     (Intercept)     I(Calc/1000) I((Calc/1000)^2) 
      0.98638196       0.22952885      -0.01899723

 

Les représentations graphiques se font en définissant les fonctions 'Rdt.Phos' et 'Rdt.Calc' (rendements en fonction respectivement de la quantité de phosphore et de calcium) puis en les traçant à l'aide de la fonction 'plot'. Les trois points moyens correspondants sont ajoutés grâce à la fonction 'points'. La fonction 'par' permet de modifier certains paramètres graphiques, en l'occurence de partitionner la fenêtre en deux (l'agrandir ensuite).

par(mfrow=c(1,2))

Rdt.Phos <- function(x) model3$coef[1] + model3$coef[2]*x + model3$coef[3]*x^2
plot(Rdt.Phos,xlim=c(0.1,0.3),ylim=c(1.1,1.8), 
xlab="Phosphore (t/ha)", ylab="Rendements (t/ha)",las=1)
points(c(0.1,0.2,0.3),
c(mean(Rend[Phos==100])/1.8,mean(Rend[Phos==200])/1.8,mean(Rend[Phos==300])/1.8), type="p")

Rdt.Calc <- function(x) model4$coef[1] + model4$coef[2]*x + model4$coef[3]*x^2
plot(Rdt.Calc,xlim=c(1,8),ylim=c(1.1,1.8), 
xlab="Calcium (t/ha)", ylab="Rendements (t/ha)",las=1)
points(c(1,4.5,8),
c(mean(Rend[Calc==1000])/1.8,mean(Rend[Calc==4500])/1.8,mean(Rend[Calc==8000])/1.8),type="p")

exp0653.jpg (61 Ko)

 

    2° Rendement maximum

Les quantités de phosphore et de calcium correspondant aux maximums des deux courbes s'obtiennent en dérivant les expressions du rendement et en annulant les dérivées obtenues. Si on refuse jusqu'au bout de prendre une feuille de papier et un crayon et que l'on ne sait même plus dériver un polynôme, cela peut se faire de la façon suivante.

Pour le phosphore :

round(model3$coef,digit=4)  # rappel des coefficients
     (Intercept)     I(Phos/1000) I((Phos/1000)^2) 
          0.2864           9.9259         -16.9753 
D(expression(0.2864+9.9259*P-16.9753*P^2),"P")   # dérivée du rendement
9.9259 - 16.9753 * (2 * P)
solve(16.9753*2,9.9259)   # détermination de la valeur qui annule la dérivée
[1] 0.292363


Pour le calcium :

round(model4$coef,digit=4)
     (Intercept)     I(Calc/1000) I((Calc/1000)^2) 
          0.9864           0.2295          -0.0190 
D(expression(0.9864+0.2295*C-0.0190*C^2),"C")
0.2295 - 0.019 * (2 * C)
solve(0.019*2,0.2295)
[1] 6.039474

 

    3° Fumure optimale

Si 1 kg d'engrais phosphorique est rentabilisé par un gain de 3 kg de blé, la quantité optimale de phosphore correspond au point où la dérivée du rendement est égale à 3. En reprenant l'expression de la dérivée du rendement par rapport au phosphore, obtenue ci-dessus, on obtient :

solve(16.9753*2, 9.9259-3)
[1] 0.2039993

De même, si 1 kg de chaux est rentabilisé par un gain de 0,1 kg de blé, la quantité optimale de chaux vaut :

solve(0.019*2, 0.2295-0.1)
[1] 3.407895

 

    4° Deux facteurs considérés simultanément

Il faut réaliser ici une régression multiple faisant intervenir les deux facteurs, leurs carrés et leur produit qui traduit l'interaction :

Rdt <- Rend/1.8 ; P <- Phos/1000 ; Ca <- Calc/1000 
model5 <- lm(Rdt ~ P + Ca + I(P^2) + I(Ca^2) + P*Ca)
round(model5$coef, digit=5)
(Intercept)           P          Ca      I(P^2)     I(Ca^2)        P:Ca 
   -0.27805    10.28307     0.24540   -16.97531    -0.01900    -0.07937

L'équation modifiée (rendement diminué des coûts) s'obtient en diminuant le coefficient de P de 3 unités et en diminuant celui de Ca de 0,1.

Il est alors possible de déterminer le rendement maximum (respectivement la fumure optimale) en annulant les dérivées partielles du rendement (respectivement du rendement diminué des coûts) par rapport aux deux facteurs et en résolvant le système obtenu.

Voir les calculs

 

    5° Courbes d'isoréponse

Les courbes d'isoréponse projetées sur le plan (P, Ca) s'obtiennent  à l'aide de la fonction 'contour'. Il faut préalablement définir une série de points pour le phosphore, une série de points pour le calcium et enfin une matrice contenant les rendements pour tous les couples (phosphore, calcium) définis par les deux séries précédentes. On peut rajouter des options, comme par exemple un dégradé de couleurs en fonction du rendement, à l'aide de la fonction 'image'.

Pour le rendement brut :

x <- seq(0.1,0.3,length=100)
y <- seq(1,8,length=100)

RdtBrut <- function(x,y) -0.2788+(10.29*x)+(0.2454*y)-(16.99*x^2)-(0.019*y^2)-(0.07929*x*y)
z <- outer(x,y,RdtBrut)

image(x,y,z,col=heat.colors(500),xlab="Phosphore (t/ha)", ylab="Calcium (t/ha)", las=1)
box()
contour(x,y,z,labcex=0.9,levels=c(1.2,1.4,1.6,1.8,1.9),add=T)

exp0654.jpg (67 Ko)

Pour le rendement diminué des coûts de la fumure :

La procédure est exactement la même sauf que l'on travaille avec l'équation modifiée. On peut reprendre le même quadrillage que précédemment mais il faut redéfinir la matrice des rendements diminués des coûts ainsi que les niveaux souhaités :

RdtModif <- function(x,y) -0.2788+(7.29*x)+(0.1454*y)-(16.99*x^2)-(0.019*y^2)-(0.07929*x*y)
zz <- outer(x,y,RdtModif)

image(x,y,zz,col=heat.colors(500),xlab="Phosphore (t/ha)", ylab="Calcium (t/ha)", las=1)
box()
contour(x,y,zz,labcex=0.9,levels=c(0.4,0.6,0.7),add=T)

wpe3.jpg (18354 octets)

 

Notons qu'il est possible de représenter la surface de réponse qui correspond à l'équation de régression, dans l'espace à trois dimensions (P, Ca, Rdt), et de représenter les différents niveaux de rendement. Une telle représentation dans l'espace est cependant peu exploitable.

Haut de la page
Autres illustrations avec R
Voir la représentation



Dernière mise à jour : janvier 2005