TP2-3
Modélisation de la taille du cerveau des chauves-souris

Chargement des packages

library(tidyverse)
library(corrplot)
library(emmeans)
library(ppcor)
library(car)
theme_set(theme_minimal())

Étude du lien entre le poids du cerveau et la masse corporelle

1

df_bats <- 
  read.table("bats.csv", header = TRUE, sep = ";", skip = 3) %>%
  mutate(Diet = factor(Diet))

df_phyto <- filter(df_bats, Diet == "1")

2

ggplot(df_phyto) +
   aes(x = BOW, y = BRW) +
   geom_point() +
   labs(x = "Masse totale du cerveau", y = "Masse corporelle")

3

On cherche à modéliser la masse du cerveau en fonction de la masse corporelle. Le modèle mathématique associé est donc

\[Y_i = a + b x_i + E_i,\]

où

• \(Y_i\) est la masse du cerveau de l’individu \(i\),
• \(x_i\) est la masse corporelle de l’individu \(i\),
• \(\left\{E_i\right\} \overset{iid}{\sim} \mathcal{N}\left(0, \sigma^2\right)\),
• \(1 \leq i \leq n = 29\).

4

reg_simple <- lm(BRW ~ BOW, data = df_phyto)
reg_simple

## 
## Call:
## lm(formula = BRW ~ BOW, data = df_phyto)
## 
## Coefficients:
## (Intercept)          BOW  
##       623.4          9.0

par(mfrow = c(2, 2)) 
plot(reg_simple)

5

On observe que le point 7 est atypique avec une distance de Cook supérieure à 1. En vérifiant avec cooks.distance(reg_simple), on s’aperçoit qu’elle vaut 100.08. On décide donc de la retirer des observations.

df_phyto2 <- df_phyto[-7, ]
reg_simple2 <- lm(BRW ~ BOW, data = df_phyto2)
par(mfrow = c(2, 2)) 
plot(reg_simple2)

Les graphes sont nettement plus proches de ce que l’on peut attendre lorsqu’on a des hypothèses respectées. Notons toutefois qu’on observe une tendance dans la moyenne (premier graphe) et une tendance dans la variance (troisième graphe) pour les petites valeurs ajustées. La normalité est validée et le point \(6\) est atypique mais nous allons le garder pour la suite de l’analyse.

Dans une analyse réelle qui sort du cadre formel de ce TP, on retirerait itérativement les points atypiques, à commencer par le point \(6\).

6

summary(reg_simple2)

## 
## Call:
## lm(formula = BRW ~ BOW, data = df_phyto2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -269.76  -93.33    8.73  112.93  322.55 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 346.5452    35.4920   9.764 3.48e-10 ***
## BOW          14.5099     0.4285  33.860  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 141.8 on 26 degrees of freedom
## Multiple R-squared:  0.9778, Adjusted R-squared:  0.977 
## F-statistic:  1147 on 1 and 26 DF,  p-value: < 2.2e-16

On effectue le test \(H_0 = \{b = 0\}\) contre \(H_1 = \{b \neq 0\}\). La statistique de test est \(F_{obs}= 1146.5\) et on a \(p<0.05\). On peut donc rejeter l’hypothèse nulle : notre modèle est significatif.

L’ordonnée à l’origine vaut \(346.5452\) et la pente \(14.5099>0\) donc plus la masse corporelle augmente, plus la masse totale du cerveau augmente.

On note que le coefficient de determination \(R^2\) vaut \(0.9778\) donc le modèle a un fort pouvoir prédictif.

Étude de la contribution des différentes parties du cerveau à sa masse totale

1

df_phyto_num <- select_if(df_phyto2, is.numeric)
df_phyto_num %>% 
  cor() %>% 
  corrplot()

Toutes les variables sont positivement corrélées entre elles.

df_phyto_num %>% 
  pcor() %>% 
  getElement("estimate") %>% 
  # `$`(estimate) %>% # alternative à getElement()
  corrplot()

Cependant, si on s’intéresse à la corrélation partielle (quand on prend en compte l’effet des autres variables), le sens de la corrélation s’inverse entre plusieurs couples.

2

On cherche à modéliser la masse du cerveau en fonction d’autres variables quantitatives. Le modèle mathématique associé est donc

\[Y_i = \beta_0 + \beta_1 x_{i,1} + \beta_2 x_{i,2} + \beta_3 x_{i,3} + E_i,\]

où

• \(Y_i\) est la masse du cerveau de l’individu \(i\),
• \(x_{i,1}\) est le volume de la partie auditive de l’individu \(i\),
• \(x_{i,2}\) est le volume de la zone olfactive de l’individu \(i\),
• \(x_{i,3}\) est le volume de l’hippocampe de l’individu \(i\),
• \(\left\{E_i\right\} \overset{iid}{\sim} \mathcal{N}\left(0, \sigma^2\right)\),
• \(1 \leq i \leq n = 29\).

reg_multiple <- lm(BRW ~ AUD + MOB + HIP, data = df_phyto2)

par(mfrow = c(2, 2)) 
plot(reg_multiple)

• On observe une tendance dans la moyenne.
• On valide la normalité.
  • De plus, on ne peut pas rejeter l’hypothèse nulle du test de Shapiro shapiro.test(reg_multiple$residuals).
• On observe une tendance dans la variance.
• Le point 6 est atypique.

3

reg_0 <- lm(BRW ~ 1, data = df_phyto2)
anova(reg_0, reg_multiple)

## Analysis of Variance Table
## 
## Model 1: BRW ~ 1
## Model 2: BRW ~ AUD + MOB + HIP
##   Res.Df      RSS Df Sum of Sq      F    Pr(>F)    
## 1     27 23562448                                  
## 2     24   603265  3  22959183 304.47 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

On effectue le test \(H_0 = \{\forall i\geq1, \beta_i = 0\}\) contre \(H_1 = \{\exists i\geq1, \beta_i \neq 0\}\). La statistique de test est \(F_{obs}= 304.47\) et on a \(p<0.05\). On peut donc rejeter l’hypothèse nulle : notre modèle est significatif et au moins une des variables a de l’influence.

On note par ailleurs que

• la somme des carrés totale vaut 23 562 448 sur 27 degrés de liberté,
• la somme des carrés résiduelle vaut 603 265 sur 24 degrés de liberté,
• la somme des carrés due au modèle vaut 22 959 183 sur 3 degrés de liberté.

4

anova(reg_multiple) # Type I

## Analysis of Variance Table
## 
## Response: BRW
##           Df   Sum Sq  Mean Sq F value    Pr(>F)    
## AUD        1  6817133  6817133 271.210 1.397e-14 ***
## MOB        1 15409397 15409397 613.040 < 2.2e-16 ***
## HIP        1   732653   732653  29.148 1.519e-05 ***
## Residuals 24   603265    25136                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• AUD apporte de l’information de manière significative sur BRW,
• MOB apporte de l’information de manière significative en plus de AUD sur BRW,
• HIP apporte de l’information de manière significative en plus de AUD et MOB sur BRW.

Anova(reg_multiple) # Type II

## Anova Table (Type II tests)
## 
## Response: BRW
##            Sum Sq Df F value    Pr(>F)    
## AUD       1572875  1  62.575 3.850e-08 ***
## MOB         14152  1   0.563    0.4603    
## HIP        732653  1  29.148 1.519e-05 ***
## Residuals  603265 24                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• l’effet marginal de AUD (en plus de MOB et HIP) sur BRW est significatif,
• l’effet marginal de MOB (en plus de AUD et HIP) sur BRW n’est pas significatif,
• l’effet marginal de HIP (en plus de AUD et MOB) sur BRW est significatif.

summary(reg_multiple)

## 
## Call:
## lm(formula = BRW ~ AUD + MOB + HIP, data = df_phyto2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -268.55  -68.84    9.88   61.66  375.34 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -312.692     76.628  -4.081  0.00043 ***
## AUD           47.989      6.067   7.910 3.85e-08 ***
## MOB           -2.444      3.257  -0.750  0.46034    
## HIP           15.981      2.960   5.399 1.52e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 158.5 on 24 degrees of freedom
## Multiple R-squared:  0.9744, Adjusted R-squared:  0.9712 
## F-statistic: 304.5 on 3 and 24 DF,  p-value: < 2.2e-16

Les p-valeurs de cette table sont les mêmes que celles de l’instruction Anova() (type II).

Le coeffidient associé à AUD est \(47.989>0\), celui de MOB est \(-2.444<0\) et celui de HIP est \(15.981>0\).

MOB est négativement associé à BRW, ce qui est logique car le coefficient de corrélation partiel obtenu plus tôt l’était aussi.

5

Selection forward

step(reg_0, scope = BRW ~ AUD + MOB + HIP, direction = "forward")

## Start:  AIC=384
## BRW ~ 1
## 
##        Df Sum of Sq      RSS    AIC
## + HIP   1  21373972  2188475 319.46
## + MOB   1  20539668  3022780 328.51
## + AUD   1   6817133 16745314 376.44
## <none>              23562448 384.00
## 
## Step:  AIC=319.46
## BRW ~ HIP
## 
##        Df Sum of Sq     RSS    AIC
## + AUD   1   1571059  617417 286.03
## <none>              2188475 319.46
## + MOB   1     12335 2176140 321.30
## 
## Step:  AIC=286.03
## BRW ~ HIP + AUD
## 
##        Df Sum of Sq    RSS    AIC
## <none>              617417 286.03
## + MOB   1     14152 603265 287.38

## 
## Call:
## lm(formula = BRW ~ HIP + AUD, data = df_phyto2)
## 
## Coefficients:
## (Intercept)          HIP          AUD  
##     -277.68        13.80        47.96

Le modèle selectionné est le modèle sans la vafiable MOB.

Selection stepwise

step(reg_0, scope = BRW ~ AUD + MOB + HIP, direction = "both")

## Start:  AIC=384
## BRW ~ 1
## 
##        Df Sum of Sq      RSS    AIC
## + HIP   1  21373972  2188475 319.46
## + MOB   1  20539668  3022780 328.51
## + AUD   1   6817133 16745314 376.44
## <none>              23562448 384.00
## 
## Step:  AIC=319.46
## BRW ~ HIP
## 
##        Df Sum of Sq      RSS    AIC
## + AUD   1   1571059   617417 286.03
## <none>               2188475 319.46
## + MOB   1     12335  2176140 321.30
## - HIP   1  21373972 23562448 384.00
## 
## Step:  AIC=286.03
## BRW ~ HIP + AUD
## 
##        Df Sum of Sq      RSS    AIC
## <none>                617417 286.03
## + MOB   1     14152   603265 287.38
## - AUD   1   1571059  2188475 319.46
## - HIP   1  16127898 16745314 376.44

## 
## Call:
## lm(formula = BRW ~ HIP + AUD, data = df_phyto2)
## 
## Coefficients:
## (Intercept)          HIP          AUD  
##     -277.68        13.80        47.96

On sélectionne le même modèle.

Étude du lien entre le volume de la partie auditive et le régime alimentaire

1

ggplot(df_bats) +
  aes(x = Diet, y = AUD, fill = Diet) +
  geom_violin() +
  geom_boxplot(alpha = 0) +
  geom_point() +
  theme_minimal() +
  theme(legend.position = "none")

2

On cherche à modéliser le volume de la partie auditive en fonction régime alimentaire, une variable qualitative. Le modèle mathématique associé est donc

\[Y_{i,k} = \mu + \alpha_i + E_{i,k},\] où

• \(Y_{i,k}\) est le volume de la partie auditive du \(k^{\text{ème}}\) individu de la \(i^{\text{ème}}\) diète,
• \(\mu\) est l’effet moyen,
• \(\alpha_{i}\) est l’effet spécifique de la diète \(i\),
• \(\left\{E_{i,k}\right\} \overset{iid}{\sim} \mathcal{N}\left(0, \sigma^2\right)\),
• \(1 \leq i \leq I = 4\),
• \(1 \leq k \leq n_i\),
• \(n_1 = 29\), \(n_2 = 5\), \(n_3 = 27\), \(n_4 = 2\).

3

reg_anova <- lm(AUD ~ Diet, data = df_bats)
par(mfrow = c(2, 2)) 
plot(reg_anova)

On valide toutes les hypothèses du modèle.

4

reg_0 <- lm(AUD ~ 1, data = df_bats)
anova(reg_0, reg_anova)

## Analysis of Variance Table
## 
## Model 1: AUD ~ 1
## Model 2: AUD ~ Diet
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     62 1464.3                           
## 2     59 1398.3  3    66.069 0.9293 0.4323

On effectue le test \(H_0 = \{\alpha_1 = \alpha_2 = \alpha_3 = \alpha_4 = 0\}\) contre \(H_1 = \{\exists i, \alpha_i \neq 0\}\). La statistique de test est \(F_{obs}= 0.9293\) et on a \(p>0.05\). On ne peut donc pas rejeter l’hypothèse nulle : notre modèle n’est pas significatif.

5

On peut effectuer une analyse de la variance à un facteur pour modéliser le volume de la partie auditive en fonction de la variable Clade.

summary(lm(AUD ~ Clade, data = df_bats))

## 
## Call:
## lm(formula = AUD ~ Clade, data = df_bats)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.988 -3.003 -1.151  2.048 19.239 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.8281     1.0452   7.489 4.01e-10 ***
## CladeII       0.5328     1.4457   0.369    0.714    
## CladeIII      3.7619     4.9026   0.767    0.446    
## CladeIV      -2.3725     1.5385  -1.542    0.128    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.79 on 59 degrees of freedom
## Multiple R-squared:  0.07559,    Adjusted R-squared:  0.02859 
## F-statistic: 1.608 on 3 and 59 DF,  p-value: 0.1971

La p-valeur associée à cette modélisation est \(p=0.1971>0.05\) donc le modèle n’est pas significatif.

Étude du lien entre le volume de la partie auditive, le régime alimentaire et le clade

1

table(df_bats$Clade, df_bats$Diet)

##      
##        1  2  3  4
##   I   10  4  7  0
##   II  19  1  1  2
##   III  0  0  1  0
##   IV   0  0 18  0

Le plan d’experience comporte des 0, il n’est ni orthogonal, ni équilibré, ni complet. On ne pourra pas estimer les effets d’interactions.

2

df_bats_merged <- 
  df_bats %>% 
  filter(Diet != 4) %>% 
  mutate(Diet = fct_drop(Diet),
         Clade = fct_collapse(Clade, II = c("II", "III")))

table(df_bats_merged$Clade, df_bats_merged$Diet)

##     
##       1  2  3
##   I  10  4  7
##   II 19  1  2
##   IV  0  0 18

Il y a toujours des zéros, on ne pourra pas estimer les effets d’interaction.

3

ggplot(df_bats_merged) +
  aes(x = Clade, y = AUD, fill = Clade) +
  geom_violin() +
  geom_boxplot(alpha = 0) +
  geom_point(alpha = 0.5) +
  facet_grid(~ Diet)

df_bats_merged %>% 
  group_by(Diet, Clade) %>% 
  summarise(`Mean of AUD` = mean(AUD)) %>% 
  ggplot() +
  aes(x = Diet, y = `Mean of AUD`, color = Clade) +
  geom_line(aes(group = Clade, linetype = Clade)) +
  geom_point()

4

On cherche à modéliser le volume de la partie auditive en fonction de l’évolution et du régime aliementaire. On ne peut pas estimer les interactions entre ces variables car le dispositif est incomplet. Le modèle mathématique associé est donc

\[Y_{i,j,k} = \mu + \alpha_i + \beta_j + E_{i,j,k}\] où

• \(Y_{i,k}\) est le volume de la partie auditive du \(k^{\text{ème}}\) individu du \(i^{\text{ème}}\) clade et de la \(j^{\text{ème}}\) diète,
• \(\mu\) est l’effet moyen,
• \(\alpha_{i}\) est l’effet spécifique du clade \(i\),
• \(\beta{j}\) est l’effet spécifique de la diète \(j\),
• \(\left\{E_{i,j,k}\right\} \overset{iid}{\sim} \mathcal{N}\left(0, \sigma^2\right)\),
• \(1 \leq i \leq I = 3\),
• \(1 \leq j \leq J = 3\),
• \(1 \leq k \leq n_{i,j}\),
• \(n_{1,1} = 10\), \(n_{1,2} = 4\), \(n_{1,3} = 7\), \(n_{2,1} = 19\), \(n_{2,2} = 1\), \(n_{2,3} = 2\), \(n_{3,1} = n_{3,2} = 0\), \(n_{3,3} = 18\).

5

reg_anova2 <- lm(AUD ~ Clade + Diet, data = df_bats_merged)

par(mfrow = c(2, 2)) 
plot(reg_anova2)

On valide toutes les hypothèses.

summary(reg_anova2)

## 
## Call:
## lm(formula = AUD ~ Clade + Diet, data = df_bats_merged)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.305 -3.116 -1.031  1.883 19.126 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.9971     1.3900   5.753 3.82e-07 ***
## CladeII       0.4766     1.6198   0.294   0.7697    
## CladeIV      -3.6897     2.0061  -1.839   0.0712 .  
## Diet2        -2.8964     2.4543  -1.180   0.2429    
## Diet3         1.1481     1.9734   0.582   0.5630    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.834 on 56 degrees of freedom
## Multiple R-squared:  0.1022, Adjusted R-squared:  0.03811 
## F-statistic: 1.594 on 4 and 56 DF,  p-value: 0.1885

anova(reg_anova2)

## Analysis of Variance Table
## 
## Response: AUD
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Clade      2   95.84  47.922  2.0506 0.1382
## Diet       2   53.19  26.597  1.1381 0.3277
## Residuals 56 1308.72  23.370

Anova(reg_anova2)

## Anova Table (Type II tests)
## 
## Response: AUD
##            Sum Sq Df F value Pr(>F)
## Clade       88.46  2  1.8926 0.1602
## Diet        53.19  2  1.1381 0.3277
## Residuals 1308.72 56

Le modèle n’est pas significatif.

Volume de la partie auditive, régime alimentaire et masse du cerveau

1

ggplot(df_bats_merged) +
  aes(x = BRW, y = AUD, color = Diet, shape = Diet) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

2

On cherche à modéliser le volume de la partie auditive en fonction du régime aliementaire et du poids du cerveau total. Le modèle mathématique associé est donc

\[Y_{i,k} = \mu + \alpha_i + (\beta + \gamma_i) \times x_{i,k} + E_{i,k}\]

où

• \(Y_{i,k}\) est le volume de la partie auditive du \(k^{\text{ème}}\) individu de la \(i^{\text{ème}}\) diète,
• \(\mu\) est l’effet moyen,
• \(\alpha_{i}\) est l’effet spécifique de la diète \(i\),
• \(\beta\) est la pente commune,
• \(\gamma_{i}\) est la pente spécifique de la diète \(i\),
• \(x_{i,k}\) est le poids du cerveau total du \(k^{\text{ème}}\) individu de la \(i^{\text{ème}}\) diète,
• \(\left\{E_{i,k}\right\} \overset{iid}{\sim} \mathcal{N}\left(0, \sigma^2\right)\),
• \(1 \leq i \leq I = 3\),
• \(1 \leq k \leq n_{i}\),
• \(n_{1} = 29\), \(n_{2} = 5\), \(n_{3} = 27\).

3

reg_ancova <- lm(AUD ~ Diet * BRW, data = df_bats_merged)
par(mfrow = c(2, 2)) 
plot(reg_ancova)

Le point 7 est atypique. Dans le cadre formel de ce TP, on continue l’analyse sans le retirer en validant les hypothèses de ce modèle.

4

reg_0 <- lm(AUD ~ 1, data = df_bats_merged)
anova(reg_0, reg_ancova)

## Analysis of Variance Table
## 
## Model 1: AUD ~ 1
## Model 2: AUD ~ Diet * BRW
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     60 1457.76                                  
## 2     55  661.21  5    796.55 13.252 1.767e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Le modèle est significatif.

On note par ailleurs que

• la somme des carrés totale vaut 1457.76 sur 60 degrés de liberté,
• la somme des carrés résiduelle vaut 661.21 sur 55 degrés de liberté,
• la somme des carrés due au modèle vaut 796.55 sur 5 degrés de liberté.

anova(reg_ancova)

## Analysis of Variance Table
## 
## Response: AUD
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Diet       2  60.58   30.29  2.5195   0.08976 .  
## BRW        1 384.52  384.52 31.9844 5.756e-07 ***
## Diet:BRW   2 351.46  175.73 14.6173 8.109e-06 ***
## Residuals 55 661.21   12.02                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(reg_ancova)

## Anova Table (Type II tests)
## 
## Response: AUD
##           Sum Sq Df F value    Pr(>F)    
## Diet        9.22  2  0.3836    0.6832    
## BRW       384.52  1 31.9844 5.756e-07 ***
## Diet:BRW  351.46  2 14.6173 8.109e-06 ***
## Residuals 661.21 55                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dans les deux cas, les effets d’intercation sont significatifs, ont conserve alors le modèle complet avec interactions.

summary(reg_ancova)

## 
## Call:
## lm(formula = AUD ~ Diet * BRW, data = df_bats_merged)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7924 -1.6128 -0.4486  0.6743 17.3482 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.9815910  0.8338761   7.173 1.96e-09 ***
## Diet2       -3.9611060  3.1610203  -1.253  0.21547    
## Diet3       -4.2154266  1.3218640  -3.189  0.00236 ** 
## BRW          0.0016506  0.0003758   4.393 5.15e-05 ***
## Diet2:BRW    0.0066103  0.0068400   0.966  0.33806    
## Diet3:BRW    0.0098132  0.0018410   5.330 1.89e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.467 on 55 degrees of freedom
## Multiple R-squared:  0.5464, Adjusted R-squared:  0.5052 
## F-statistic: 13.25 on 5 and 55 DF,  p-value: 1.767e-08

5

df_bats_merged %>%
  group_by(Diet) %>%
  summarise(Effectif = n(),
            `Mean of BRW` = mean(BRW),
            `Mean of AUD` = mean(AUD))

## # A tibble: 3 x 4
##   Diet  Effectif `Mean of BRW` `Mean of AUD`
##   <fct>    <int>         <dbl>         <dbl>
## 1 1           29         1410.          8.31
## 2 2            5          384.          5.20
## 3 3           27          432.          6.72

lsmeans(reg_ancova, pairwise ~ Diet, adjust = "bonferroni")

## NOTE: Results may be misleading due to involvement in interactions

## $lsmeans
##  Diet lsmean    SE df lower.CL upper.CL
##  1      7.46 0.673 55     6.11      8.8
##  2      9.40 3.805 55     1.77     17.0
##  3     12.01 1.066 55     9.87     14.1
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE df t.ratio p.value
##  1 - 2       -1.94 3.86 55 -0.503  1.0000 
##  1 - 3       -4.55 1.26 55 -3.611  0.0020 
##  2 - 3       -2.61 3.95 55 -0.660  1.0000 
## 
## P value adjustment: bonferroni method for 3 tests

La diète 3 a une influence significativement différente de la diète 1, après correction pour la multiplicité par la méthode de Bonferroni. En revanche, les diètes 1 et 2 ne sont pas significativement différentes, ni les diètes et 2 et 3.

---