1 . Introduction

Vous avez étudié la hauteur de plants de maïs chez différentes variétés(graines) et avec trois types d’engrais différents. Les engrais et les variétés peuvent être appliqués uniquement à une parcelle d’un coup. Vous souhaitez à présent analyser vos résultats.

1. Vous ne vous souvenez pas très bien de votre plan d’expérience. Pour vous rappeler de celui-ci, chargez les données récoltées (TP5_data.csv) et affichez-les (\(y\equiv\) hauteur des plants de maïs [cm]). Essayez de représenter le schéma de l’expérience sur une feuille à part.

mydata <- read.csv("TP5_data.csv")
factors <- c("graine", "engrais", "parcelle", "plant", "mesure")
mydata[factors] <- lapply(mydata[factors], factor)

pander(summary(mydata))
Table continues below
y Ymoyen graine engrais parcelle plant mesure
Min. :10.09 Min. :10.36 1:96 1:96 1:72 1:72 1:144
1st Qu.:12.68 1st Qu.:12.73 2:96 2:96 2:72 2:72 2:144
Median :13.74 Median :13.67 3:96 3:96 3:72 3:72
Mean :14.00 Mean :14.00 4:72 4:72
3rd Qu.:15.63 3rd Qu.:15.62
Max. :17.52 Max. :17.07
NA’s :252
IDparcelle
Min. :111.0
1st Qu.:131.8
Median :222.5
Mean :222.5
3rd Qu.:313.2
Max. :334.0 
ggplot(mydata,aes(x=engrais,y=y)) + geom_point(alpha = 0.5, size=3) + facet_wrap(~graine) +
               labs(x= "Engrais",
                    y= "Hauteur du plant")

2. Vos données sont-elles balancées?
Oui

pander(table(mydata[,c("graine", "engrais","parcelle","plant")]))
plant 1 2 3 4
graine engrais parcelle
1 1 1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
2 1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
3 1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
2 1 1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
2 1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
3 1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
3 1 1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
2 1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
3 1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2

Lors de ce TP, les modèles seront codés en utilisant la méthode ‘sum contrast’

Avec la méthode ‘treatment/dummy contrast’, les paramètres estimés sont \(\mu + \alpha_1\), \(\alpha_2 - \alpha_1\) et \(\alpha_3 - \alpha_1\). Cette méthode est surtout utile lorsque l’on a un niveau du facteur qui correspond au témoin de l’expérience.
Avec la méthode ‘sum contrast’, les paramètres estimés sont \(\mu\), \(\alpha_1\) et \(\alpha_2\). Étant donné que \(\sum \hat{\alpha}_i = 0\), on peut trouver facilement \(\hat{\alpha}_3=-\hat{\alpha}_1 - \hat{\alpha}_2\) si les données sont balancées.

2 . Modèle 1

3. Pourquoi est-ce intéressant d’inclure les facteurs parcelle et plant dans l’analyse des résultats ?

Pour éviter que la part de variance de la réponse due à ces facteurs ne soit confondue avec la variance des résidus. Si c’est le cas, on observera moins facilement l’effet des facteurs \(\alpha\) et \(\beta\).

4. Combien de réplications ont été faites pour chaque niveau de facteur ainsi que pour la combinaison des deux facteurs (traitement)? Les deux mesures sur un même plant sont-elles considérées comme des réplications?

Rappel:
- Unité expérimentale: la plus petite unité à laquelle un traitement est appliqué
- Réplication: consiste à appliquer chaque traitement à plusieurs unités expérimentales et mutuellement indépendantes
- Pseudo-réplication: consiste à appliquer plusieurs mesures dans une même unité expérimentale

12 réplications pour le facteur ‘graine’ et pour le facteur ‘engrais’ (car ils sont appliqués par parcelle). Pour la combinaison des deux facteurs, nous avons 4 réplications.

Les 4 plantes ainsi que les deux mesures par plante sont considérées comme des pseudo-réplications.

table(mydata$engrais)/8 # 4 plantes par parcelle et deux mesures, donc 8 pseudoréplications 

 1  2  3 
12 12 12 
table(mydata[,c("graine", "engrais")])/8 # interaction
      engrais
graine 1 2 3
     1 4 4 4
     2 4 4 4
     3 4 4 4

5. Faites un modèle en incluant les facteurs graine, engrais, parcelle et plant de manière judicieuse. Affichez le résumé de ce modèle et faites un test ANOVA. Ecrivez l’équation de votre modèle.

(Pour avoir les pvaleurs des effets fixes et l’anova correcte, il faut charger le package lmerTest)

\(y_{ijklm} = \mu + \alpha_i + \beta_j + \alpha\beta_{ij} + c_{ijk} + d_{l(ijk)} + \epsilon_{ijklm}\)

mod1 <- lmer(y~graine*engrais + (1|IDparcelle) + (1|IDparcelle:plant), data = mydata, contrasts = list(graine=c("contr.sum"),engrais=c("contr.sum")))
summary(mod1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: y ~ graine * engrais + (1 | IDparcelle) + (1 | IDparcelle:plant)
   Data: mydata

REML criterion at convergence: -183.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.50626 -0.48721  0.01317  0.48128  2.30535 

Random effects:
 Groups           Name        Variance Std.Dev.
 IDparcelle:plant (Intercept) 0.072281 0.26885 
 IDparcelle       (Intercept) 2.184798 1.47811 
 Residual                     0.002067 0.04547 
Number of obs: 288, groups:  IDparcelle:plant, 144; IDparcelle, 36

Fixed effects:
                  Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)      14.000013   0.247382 26.999978  56.593  < 2e-16 ***
graine1          -0.910799   0.349851 26.999974  -2.603  0.01482 *  
graine2           0.006897   0.349851 26.999974   0.020  0.98442    
engrais1         -1.083225   0.349851 26.999974  -3.096  0.00453 ** 
engrais2         -0.029421   0.349851 26.999974  -0.084  0.93360    
graine1:engrais1  0.420422   0.494764 27.000006   0.850  0.40294    
graine2:engrais1 -0.767448   0.494764 27.000006  -1.551  0.13251    
graine1:engrais2  0.223478   0.494764 27.000006   0.452  0.65510    
graine2:engrais2 -0.195061   0.494764 27.000006  -0.394  0.69649    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) grain1 grain2 engrs1 engrs2 grn1:1 grn2:1 grn1:2
graine1      0.000                                                 
graine2      0.000 -0.500                                          
engrais1     0.000  0.000  0.000                                   
engrais2     0.000  0.000  0.000 -0.500                            
gran1:ngrs1  0.000  0.000  0.000  0.000  0.000                     
gran2:ngrs1  0.000  0.000  0.000  0.000  0.000 -0.500              
gran1:ngrs2  0.000  0.000  0.000  0.000  0.000 -0.500  0.250       
gran2:ngrs2  0.000  0.000  0.000  0.000  0.000  0.250 -0.500 -0.500
anova(mod1)
Type III Analysis of Variance Table with Satterthwaite's method
                  Sum Sq   Mean Sq NumDF DenDF F value   Pr(>F)   
graine         0.0185404 0.0092702     2    27  4.4845 0.020816 * 
engrais        0.0271605 0.0135803     2    27  6.5695 0.004735 **
graine:engrais 0.0090733 0.0022683     4    27  1.0973 0.377948   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6. À quoi correspondent les paramètres que vous obtenez dans la colonne estimate des effets fixes?
L’estimation de la moyenne de la réponse \(\mu\) ainsi que les estimations des \(\alpha_i\), \(\beta_j\) et \(\alpha \beta_{ij}\) pour les deux premiers niveaux de chaque facteur.

7. Combien de paramètres semblent significatifs dans votre modèle (pour un \(\alpha\) de 0.05) ?
Trois. Si on regarde l’anova, on remarque aussi que les facteurs ‘graine’ et ‘engrais’ sont significatifs, mais pas l’interaction.

8. Que vaut la covariance entre les hauteurs de deux plants de maïs d’une même parcelle?
2.184798

9. Calculez et affichez la matrive V. Pourquoi ne retrouve-t-on que deux valeurs différentes dans cette matrice? À quoi correspondent ces valeurs?
On pourrait croire qu’il n’y a que deux valeurs (mis à part les zéros) car la variance des mesures est très faible.
2.259146 = \(\sigma^2_Y = \sigma^2_c + \sigma^2_d + \sigma^2_{\epsilon}\)
2.257078 = \(\sigma^2_c + \sigma^2_d\)
2.184798 = \(\sigma^2_c\)

G <- as.matrix(getME(mod1, "Lambda") * sigma(mod1)) ** 2
Z <- as.matrix(getME(mod1,"Z"))
R <-  sigma(mod1)**2 * diag(nrow(mydata))
V <- Z %*% G %*% t(Z) + R
mybreaks = sort(unique(c(V))+0.01)
plot(V[1:16,1:16], border=NA, breaks=mybreaks)

3 . Modèle 2

10. La colonne Ymoyen correspond à la moyenne des hauteurs par parcelle pour chaque traitement. Estimez un deuxième modèle en utilisant Ymoyen comme réponse. Affichez le résumé de ce modèle.Astuce: pour utiliser la colonne Ymoyen en retirant les cases vides, vous pouvez utiliser : mydata[!is.na(mydata$Ymoyen),]

mod2 <- lm(Ymoyen~graine*engrais, data = mydata[!is.na(mydata$Ymoyen),], contrasts = list(graine=c("contr.sum"),engrais=c("contr.sum")))
summary(mod2)

Call:
lm(formula = Ymoyen ~ graine * engrais, data = mydata[!is.na(mydata$Ymoyen), 
    ], contrasts = list(graine = c("contr.sum"), engrais = c("contr.sum")))

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3713 -0.7061  0.0591  0.6658  2.8002 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      14.000013   0.247382  56.593  < 2e-16 ***
graine1          -0.910799   0.349851  -2.603  0.01482 *  
graine2           0.006897   0.349851   0.020  0.98442    
engrais1         -1.083225   0.349851  -3.096  0.00453 ** 
engrais2         -0.029421   0.349851  -0.084  0.93360    
graine1:engrais1  0.420422   0.494764   0.850  0.40294    
graine2:engrais1 -0.767448   0.494764  -1.551  0.13251    
graine1:engrais2  0.223478   0.494764   0.452  0.65510    
graine2:engrais2 -0.195061   0.494764  -0.394  0.69649    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.484 on 27 degrees of freedom
Multiple R-squared:  0.4953,    Adjusted R-squared:  0.3458 
F-statistic: 3.312 on 8 and 27 DF,  p-value: 0.009122
anova(mod2)
Analysis of Variance Table

Response: Ymoyen
               Df Sum Sq Mean Sq F value   Pr(>F)   
graine          2 19.760  9.8798  4.4845 0.020816 * 
engrais         2 28.947 14.4733  6.5695 0.004735 **
graine:engrais  4  9.670  2.4175  1.0973 0.377947   
Residuals      27 59.484  2.2031                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

11. Comparez ce modèle au modèle 1. Qu’observez-vous? Selon vous, à quelle particularité du plan d’expérience est due ce que vous observez ? Quel pourrait être l’intérêt de travailler sur les moyennes?

On obtient les mêmes résultats que pour la partie fixe du modèle 1 car le plan est balancé. C’est plus simple de travailler avec les moyennes si le plan est balancé et que l’on est pas intéressé de connaitre les composantes de variance.