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
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) |
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, |
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)) |
[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) |
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 |
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)) |
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 |
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
|
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") |
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) |
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) |
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 |
Dernière mise à jour : janvier 2005