Bases du modèle linéaire

class: center, middle, inverse, title-slide

# Bases du modèle linéaire
## Chapitre IV - ANOVA à deux facteurs
### Antoine Bichat - Émilie Lebarbier<br>AgroParisTech

---

# Chargement des données

```r
colza <- read.table("Colza.txt", header = T)
str(colza)
```

```
'data.frame':	60 obs. of  3 variables:
 $ Fertilisation: int  1 1 1 1 1 1 1 1 1 1 ...
 $ Rotation     : Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...
 $ PdsGrains    : num  27.6 16.3 11.4 38.2 38.1 24.7 22.7 21.7 20.6 19.8 ...
```

# Transformation en facteur

```r
colza$Fertilisation <- as.factor(colza$Fertilisation)
str(colza)
```

```
'data.frame':	60 obs. of  3 variables:
 $ Fertilisation: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
 $ Rotation     : Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...
 $ PdsGrains    : num  27.6 16.3 11.4 38.2 38.1 24.7 22.7 21.7 20.6 19.8 ...
```

---

# Tables

```r
table(colza$Fertilisation)
```

```

1  2 
30 30 
```

```r
table(colza$Rotation)
```

```

A  B  C 
20 20 20 
```

```r
table(colza$Fertilisation, colza$Rotation)
```

```
   
     A  B  C
  1 10 10 10
  2 10 10 10
```

---
# Description des données

```r
summary(colza$PdsGrains)
```

```
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3.40   18.00   23.45   24.02   29.10   44.10 
```

```r
sd(colza$PdsGrains)
```

```
[1] 8.936966
```

```r
by(colza$PdsGrains, colza$Fertilisation, mean)
```

```
colza$Fertilisation: 1
[1] 25.58333
-------------------------------------------------------- 
colza$Fertilisation: 2
[1] 22.46667
```

```r
by(colza$PdsGrains, colza$Rotation, mean)
```

```
colza$Rotation: A
[1] 19.96
-------------------------------------------------------- 
colza$Rotation: B
[1] 21.92
-------------------------------------------------------- 
colza$Rotation: C
[1] 30.195
```

---
# Description des données avec {tidyverse}

```r
library(tidyverse)

colza %>% 
  group_by(Fertilisation, Rotation) %>% 
  summarise(N = n(), Mean = mean(PdsGrains), 
            Min = min(PdsGrains), Median = median(PdsGrains), 
            Max = max(PdsGrains), Var = sd(PdsGrains)^2) %>% 
  arrange(Mean)
```

```
# A tibble: 6 x 8
# Groups:   Fertilisation [2]
  Fertilisation Rotation     N  Mean   Min Median   Max   Var
  <fct>         <fct>    <int> <dbl> <dbl>  <dbl> <dbl> <dbl>
1 2             A           10  15.8   3.4   16.8  27.8  55.3
2 2             B           10  19.8   8.7   16.8  34.6  68.4
3 1             B           10  24    14     22.0  39.2  54.3
4 1             A           10  24.1  11.4   22.2  38.2  74.2
5 1             C           10  28.6  20.3   28.0  40.6  34.4
6 2             C           10  31.8  23.4   31.4  44.1  52.5
```

---
# Boxplots

```r
boxplot(colza$PdsGrains ~ colza$Fertilisation * colza$Rotation,
        xlab = "Traitement", ylab = "Poids des grains")
```

---
# Graphes d'interaction

.pull-left[

```r
interaction.plot(colza$Rotation, 
                 colza$Fertilisation,
                 colza$PdsGrains,
                 fixed = TRUE, 
                 col = 2:3, 
                 leg.bty = "o")
```

<img src="chap4_files/figure-html/unnamed-chunk-5-1.png" width="504" style="display: block; margin: auto;" />
]

.pull-right[

```r
interaction.plot(colza$Fertilisation, 
                 colza$Rotation,
                 colza$PdsGrains,
                 fixed = TRUE, 
                 col = 2:4, 
                 leg.bty = "o")
```

<img src="chap4_files/figure-html/unnamed-chunk-6-1.png" width="504" style="display: block; margin: auto;" />
]

---
# Graphes avec {ggplot2}

```r
ggplot(colza) +
  aes(x = Fertilisation, y = PdsGrains, fill = Rotation) +
  geom_boxplot() +
  theme_minimal() + labs(y = "Poids des grains")
```

---

```r
colza %>% 
  group_by(Rotation, Fertilisation) %>% 
  summarise(Mean = mean(PdsGrains)) %>% 
  ggplot() +
  aes(x = Rotation, y = Mean, group = Fertilisation, 
      color = Fertilisation, linetype = Fertilisation) +
  geom_line(size = 2) +
  theme_minimal() + labs(y = "Poids moyen des grains")
```

---
# Modèle à deux facteurs avec interaction

```r
mod_2F_I <- lm(PdsGrains ~ Fertilisation * Rotation, data = colza)
mod_2F_I
```

```

Call:
lm(formula = PdsGrains ~ Fertilisation * Rotation, data = colza)

Coefficients:
             (Intercept)            Fertilisation2  
                   24.11                     -8.30  
               RotationB                 RotationC  
                   -0.11                      4.53  
Fertilisation2:RotationB  Fertilisation2:RotationC  
                    4.14                     11.41  
```

---
# Graphes de diagnostic

```r
par(mfrow = c(2, 2))
plot(mod_2F_I)
```

---
# Test du modèle

```r
mod_0 <- lm(PdsGrains ~ 1, data = colza)
anova(mod_0, mod_2F_I)
```

```
Analysis of Variance Table

Model 1: PdsGrains ~ 1
Model 2: PdsGrains ~ Fertilisation * Rotation
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1     59 4712.3                                  
2     54 3052.5  5    1659.8 5.8726 0.0002106 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---
# Test des facteurs

```r
anova(mod_2F_I)
```

```
Analysis of Variance Table

Response: PdsGrains
                       Df  Sum Sq Mean Sq F value    Pr(>F)    
Fertilisation           1  145.70  145.70  2.5776 0.1142193    
Rotation                2 1180.48  590.24 10.4417 0.0001466 ***
Fertilisation:Rotation  2  333.63  166.82  2.9511 0.0607686 .  
Residuals              54 3052.47   56.53                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

```r
car::Anova(mod_2F_I)
```

```
Anova Table (Type II tests)

Response: PdsGrains
                        Sum Sq Df F value    Pr(>F)    
Fertilisation           145.70  1  2.5776 0.1142193    
Rotation               1180.48  2 10.4417 0.0001466 ***
Fertilisation:Rotation  333.63  2  2.9511 0.0607686 .  
Residuals              3052.47 54                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---
# Coefficients

```r
summary(mod_2F_I)
```

```

Call:
lm(formula = PdsGrains ~ Fertilisation * Rotation, data = colza)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.710  -5.492  -1.125   3.752  15.200

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)                24.110      2.377  10.141 4.16e-14 ***
Fertilisation2             -8.300      3.362  -2.469   0.0168 *  
RotationB                  -0.110      3.362  -0.033   0.9740    
RotationC                   4.530      3.362   1.347   0.1835    
Fertilisation2:RotationB    4.140      4.755   0.871   0.3878    
Fertilisation2:RotationC   11.410      4.755   2.400   0.0199 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.518 on 54 degrees of freedom
Multiple R-squared:  0.3522,	Adjusted R-squared:  0.2923 
F-statistic: 5.873 on 5 and 54 DF,  p-value: 0.0002106
```

---
# Modèle à deux facteurs sans interaction

```r
mod_2F <- lm(PdsGrains ~ Fertilisation + Rotation, data = colza)
mod_2F
```

```

Call:
lm(formula = PdsGrains ~ Fertilisation + Rotation, data = colza)

Coefficients:
   (Intercept)  Fertilisation2       RotationB       RotationC  
        21.518          -3.117           1.960          10.235  
```

---
# Graphes de diagnostic

```r
par(mfrow = c(2, 2))
plot(mod_2F)
```

---
# Test des facteurs

```r
anova(mod_2F)
```

```
Analysis of Variance Table

Response: PdsGrains
              Df Sum Sq Mean Sq F value    Pr(>F)    
Fertilisation  1  145.7  145.70  2.4097 0.1262207    
Rotation       2 1180.5  590.24  9.7615 0.0002307 ***
Residuals     56 3386.1   60.47                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

```r
car::Anova(mod_2F)
```

```
Anova Table (Type II tests)

Response: PdsGrains
              Sum Sq Df F value    Pr(>F)    
Fertilisation  145.7  1  2.4097 0.1262207    
Rotation      1180.5  2  9.7615 0.0002307 ***
Residuals     3386.1 56                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---
# Coefficients

```r
summary(mod_2F)
```

```

Call:
lm(formula = PdsGrains ~ Fertilisation + Rotation, data = colza)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.002  -5.074  -1.428   5.287  16.682

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      21.518      2.008  10.718 3.49e-15 ***
Fertilisation2   -3.117      2.008  -1.552  0.12622    
RotationB         1.960      2.459   0.797  0.42877    
RotationC        10.235      2.459   4.162  0.00011 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.776 on 56 degrees of freedom
Multiple R-squared:  0.2814,	Adjusted R-squared:  0.2429 
F-statistic: 7.311 on 3 and 56 DF,  p-value: 0.0003203
```

---
# Comparaison des moyennes

```r
lsmeans::lsmeans(mod_2F, pairwise ~ Rotation, adjust = "bonferroni")
```

```
$lsmeans
 Rotation lsmean   SE df lower.CL upper.CL
 A          20.0 1.74 56     16.5     23.4
 B          21.9 1.74 56     18.4     25.4
 C          30.2 1.74 56     26.7     33.7

Results are averaged over the levels of: Fertilisation 
Confidence level used: 0.95

$contrasts
 contrast estimate   SE df t.ratio p.value
 A - B       -1.96 2.46 56 -0.797  1.0000 
 A - C      -10.23 2.46 56 -4.162  0.0003 
 B - C       -8.28 2.46 56 -3.365  0.0042

Results are averaged over the levels of: Fertilisation 
P value adjustment: bonferroni method for 3 tests 
```

---
class: center, middle, inverse

# Des questions ?

.footnote[Slides créées avec le package <b><a href="https://github.com/yihui/xaringan" target="_blank">xaringan</a></b>.]