Chapitre 11 Indice `B` de l’intérieur de la coquille

La démarche sera la même que celle des chapitres précédents. Il se peut qu’il y ait moins de commentaires.

Même jeu de données oeuf.csv qui contient différentes mesures dont l’ évaluation de la couleur de l’intérieur de la coquille (Yellownish index) - indice B , mesurée en 5 séances. Mêmes traitements (régimes).

La question est de savoir si les différents régimes induisent des indices de coloration jaune à l’intérieur de la coquille d’oeuf significativement différents avec le temps.

Mais chaque traitement n’ayant pas été appliqué sur tous les groupes d’oiseaux, l’ANOVA à mesures répétées ne pourrait pas être appliquée. Nous comparerons les effets des traitements séance par séance, puis à l’aide d’une figure on appréciera s’il y a une évolution de cet indice en fonction du temps.

11.1 Les données

bic <- read_csv("data/oeuf.csv")
bic <- bic %>% 
  select(seance, regime, no_oeuf, indbic) %>% 
  mutate(id = rep(1:30, 5), .before = 1) %>% 
  convert_as_factor(id, seance, regime)

Le tableau a été préalablement structuré en format long en Excel. J’ai ajouté un identifiant (id) pour les échantillons des séances.

glimpse(bic)

## Rows: 150
## Columns: 5
## $ id      <fct> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,~
## $ seance  <fct> seance 1, seance 1, seance 1, seance 1, seance 1, seance 1, se~
## $ regime  <fct> "Ba 0,25", "Ba 0,25", "Ba 0,25", "Ba 0,50", "Ba 0,50", "Ba 0,5~
## $ no_oeuf <dbl> 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3,~
## $ indbic  <dbl> 13.360, 13.840, 11.830, 13.040, 13.600, 14.880, 12.070, 12.680~

11.2 Visualisation boxplots

bxp <- ggplot(bic, aes(x = regime, y = indbic, fill = regime)) +
  geom_boxplot() +
  facet_grid(seance ~ .) +
  theme(axis.text.x = element_text(angle = 90, color = "bbick", vjust = 0.5, hjust = 1)) +
  ylab("Egg shell yellownish index") +
  theme_bw()
bxp

=> Variations notables entre les traitements pour certaines séances, surtout entre la séance 1 et les autres.

11.3 Détection des valeurs aberrantes extrêmes

bic_out <- bic %>%
  group_by(seance, regime) %>%
  identify_outliers(indbic)
bic_out

## [1] seance     regime     id         no_oeuf    indbic     is.outlier is.extreme
## <0 rows> (or 0-length row.names)

=> Pas de valeurs aberrantes extrêmes pour toutes les séances.

11.4 Conditions de l’ANOVA

11.4.1 Normalité

Si les données sont normalement distribuées, la p-value de Shapiro-Wilk doit être supérieure à 0,05 pour chaque régime.

bic %>%
  group_by(seance) %>%
  shapiro_test(indbic)

## # A tibble: 5 x 4
##   seance   variable statistic             p
##   <fct>    <chr>        <dbl>         <dbl>
## 1 seance 1 indbic       0.436 0.00000000121
## 2 seance 2 indbic       0.828 0.000221     
## 3 seance 3 indbic       0.954 0.214        
## 4 seance 4 indbic       0.961 0.320        
## 5 seance 5 indbic       0.842 0.000418

=> Cette hypothèse n’est pas respectée pour toutes les séances.

Créer des QQ-plots pour chaque point par séance

ggqqplot(bic, "indbic", facet.by = "seance")

=> Mais selon les QQ-plots seule la séance 1 présente un problème de normalité.

On explorera les données séance par séance pour palier au problème de normalité.

11.4.2 Homogénéité des variances

bic %>%
  select(seance, regime, indbic) %>% 
  group_by(seance) %>%
  levene_test(indbic ~ regime)

## # A tibble: 5 x 5
##   seance     df1   df2 statistic     p
##   <fct>    <int> <int>     <dbl> <dbl>
## 1 seance 1     9    20     1.66  0.164
## 2 seance 2     9    20     0.936 0.517
## 3 seance 3     9    20     1.18  0.357
## 4 seance 4     9    20     0.643 0.748
## 5 seance 5     9    20     1.15  0.377

=> Toutes les valeurs p sont > 0.05 => toutes les variances sont homogènes.

11.5 ANOVA à 1 facteur séance par séance

11.5.1 Séance 1

bic1 <- bic %>% filter(seance == "seance 1")

(bic1_out <- bic1 %>% 
  identify_outliers(indbic))

## # A tibble: 2 x 7
##   id    seance   regime no_oeuf indbic is.outlier is.extreme
##   <fct> <fct>    <fct>    <dbl>  <dbl> <lgl>      <lgl>     
## 1 29    seance 1 WC           2   30.0 TRUE       TRUE      
## 2 30    seance 1 WC           3   70.8 TRUE       TRUE

=> 2 observatiosn s’écartent énormémemt du nuage de point pour cette séance. Est-ce qu’il faut les exclure ?

#bic1 <- bic1 %>% filter(! id %in% c(29, 30))

11.5.1.1 Le modèle

lm1 <- lm(log(indbic) ~ regime, data = bic1)
Anova(lm1)

## Anova Table (Type II tests)
## 
## Response: log(indbic)
##           Sum Sq Df F value   Pr(>F)   
## regime    2.4579  9  4.8115 0.001667 **
## Residuals 1.1352 20                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La p-value > 0.05 => Pas de différence significative entre les effets des régimes sur cet indice à la séance 1.

shapiro_test(residuals(lm1))

## # A tibble: 1 x 3
##   variable       statistic  p.value
##   <chr>              <dbl>    <dbl>
## 1 residuals(lm1)     0.829 0.000229

=> Normalité pas respectée ! Même avec la transformation log().

=> Alternative : Kruskal-Wallis et Dunn

Test de Kurskal-Wallis

bic1 %>% kruskal_test(indbic ~ regime)

## # A tibble: 1 x 6
##   .y.        n statistic    df     p method        
## * <chr>  <int>     <dbl> <int> <dbl> <chr>         
## 1 indbic    30      20.2     9 0.017 Kruskal-Wallis

11.5.1.2 Comparaisons par paires

Test de Dunn

bic1 %>% 
  dunn_test(indbic ~ regime, p.adjust.method = "bonferroni") %>% 
  select(group1, group2, p, p.adj, p.adj.signif) %>% 
  filter(p.adj.signif != "ns")

## # A tibble: 1 x 5
##   group1  group2        p  p.adj p.adj.signif
##   <chr>   <chr>     <dbl>  <dbl> <chr>       
## 1 Ba 0,75 WC     0.000707 0.0318 *

=> Différence au niveau des effets de ces deux régimes.

Calcul des moyennes ci-dessous, pour le classement manuel (identifier les groupes). Le graphique automatique un peu complexe !

cm1 <- bic1 %>% 
  select(regime, indbic) %>% 
  group_by(regime) %>% 
  summarise_all(list(bic_moyen = mean, sd = sd)) %>% 
  data.table::setorder(-bic_moyen)
cm1

## # A tibble: 10 x 3
##    regime  bic_moyen     sd
##  * <fct>       <dbl>  <dbl>
##  1 WC           40.2 27.0  
##  2 YC           18.7  3.60 
##  3 Ba 2,5       16.1  1.08 
##  4 Ba 7,5       15.3  3.18 
##  5 Ba 5         15.2  0.348
##  6 Ba 1         14.1  3.22 
##  7 Ba 0,50      13.8  0.943
##  8 Ba 10        13.1  0    
##  9 Ba 0,25      13.0  1.05 
## 10 Ba 0,75      12.4  0.305

J’ajoute les lettres à la main :

cm1 <- cm1 %>% 
  mutate(groups = c("a", "ab", "ab", "ab", "ab", "ab", "ab", "ab", "ab", "b"))
cm1

## # A tibble: 10 x 4
##    regime  bic_moyen     sd groups
##  * <fct>       <dbl>  <dbl> <chr> 
##  1 WC           40.2 27.0   a     
##  2 YC           18.7  3.60  ab    
##  3 Ba 2,5       16.1  1.08  ab    
##  4 Ba 7,5       15.3  3.18  ab    
##  5 Ba 5         15.2  0.348 ab    
##  6 Ba 1         14.1  3.22  ab    
##  7 Ba 0,50      13.8  0.943 ab    
##  8 Ba 10        13.1  0     ab    
##  9 Ba 0,25      13.0  1.05  ab    
## 10 Ba 0,75      12.4  0.305 b

ggplot(data = cm1, mapping = aes(x = regime, y = bic_moyen)) +
  geom_bar(stat = "identity", color = "blue", fill = "grey", width = 0.6) +
  #ylim(0, 42) +
  geom_text(aes(label = groups), vjust = -0.5, size = 4) +
  xlab("Régimes") + ylab("Egg shell yellownish index") +
  theme(axis.text.x = element_text(angle = 45, color = "black", vjust = 1, hjust = 1)) +
  theme_bw()

11.5.2 Séance 2

11.5.2.1 Le modèle

bic2 <- bic %>% filter(seance == "seance 2")

(bic2_out <- bic2 %>% 
  identify_outliers(indbic))

## # A tibble: 2 x 7
##   id    seance   regime  no_oeuf indbic is.outlier is.extreme
##   <fct> <fct>    <fct>     <dbl>  <dbl> <lgl>      <lgl>     
## 1 5     seance 2 Ba 0,50       2   5.71 TRUE       FALSE     
## 2 29    seance 2 WC            2   7.92 TRUE       TRUE

=> 1 observation supossée extrême

bic2 <- bic2 %>% filter(id != 29)

lm2 <- lm(indbic ~ regime, data = bic2)
Anova(lm2)

## Anova Table (Type II tests)
## 
## Response: indbic
##            Sum Sq Df F value   Pr(>F)   
## regime    18.4897  9  4.3486 0.003392 **
## Residuals  8.9761 19                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La p-value est > 0.05 => Pas de différence significative entre les effetes des régimes.

shapiro_test(residuals(lm2))

## # A tibble: 1 x 3
##   variable       statistic p.value
##   <chr>              <dbl>   <dbl>
## 1 residuals(lm2)     0.950   0.187

=> Normalité pas respectée même avec la transformation log(). Groupes pas distingués avec la méthode alternative non paramétrique. La suppresion de l’observation extrême a amélioré la distribution.

La normalité étant respectée, on continue sans la valeur extrême.

11.5.2.2 Comparaisons par paires

cm2 <- (SNK.test(lm2, "regime", group = TRUE))$groups %>% 
  mutate(regime = rownames(.)) %>% 
  select(regime, indbic, groups) %>% 
  as_tibble()
cm2

## # A tibble: 10 x 3
##    regime  indbic groups
##    <chr>    <dbl> <chr> 
##  1 Ba 0,50   4.25 a     
##  2 Ba 7,5    3.54 ab    
##  3 Ba 5      3.01 ab    
##  4 YC        3.01 ab    
##  5 Ba 0,75   2.14 b     
##  6 Ba 2,5    2.13 b     
##  7 Ba 10     2.06 b     
##  8 Ba 0,25   1.94 b     
##  9 Ba 1      1.90 b     
## 10 WC        1.67 b

Visualisation des groupes

ggplot(data = cm2, mapping = aes(x = regime, y = indbic)) +
  geom_bar(stat = "identity", color = "blue", fill = "grey", width = 0.6) +
  #ylim(0, 5) +
  geom_text(aes(label = groups), vjust = -0.5, size = 4) +
  xlab("Régimes") + ylab("Egg shell yellownish index") +
  theme(axis.text.x = element_text(angle = 45, color = "bbick", vjust = 1, hjust = 1)) +
  theme_bw()

11.5.3 Séance 3

bic3 <- bic %>% filter(seance == "seance 3")

(bic3_out <- bic3 %>% 
    identify_outliers(indbic))

## # A tibble: 3 x 7
##   id    seance   regime  no_oeuf indbic is.outlier is.extreme
##   <fct> <fct>    <fct>     <dbl>  <dbl> <lgl>      <lgl>     
## 1 5     seance 3 Ba 0,50       2  10.4  TRUE       FALSE     
## 2 11    seance 3 Ba 1          2   9.15 TRUE       FALSE     
## 3 17    seance 3 Ba 5          2   9.47 TRUE       FALSE

=> Pas d’observations aberrantes extrêmes.

#bic3 <- bic3 %>% filter(id != ...)

11.5.3.1 Le modèle

lm3 <- lm(log(indbic) ~ regime, data = bic3)
Anova(lm3)

## Anova Table (Type II tests)
## 
## Response: log(indbic)
##           Sum Sq Df F value    Pr(>F)    
## regime    9.7056  9  9.8062 1.334e-05 ***
## Residuals 2.1994 20                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La p-value < 0.01 => Différence très significative entre les effetes d’au moins 2 régimes.

shapiro_test(residuals(lm3))

## # A tibble: 1 x 3
##   variable       statistic p.value
##   <chr>              <dbl>   <dbl>
## 1 residuals(lm3)     0.944   0.115

=> Normalité pas satisfaite. La transformation log() a résolu le problème.

11.5.3.2 Comparaisons par paires

cm3 <- (SNK.test(lm3, "regime", group = TRUE))$groups %>% 
  mutate(regime = rownames(.)) 
colnames(cm3)[1] <- "indbic"
cm3

##             indbic groups  regime
## Ba 5    2.00733090      a    Ba 5
## Ba 0,50 1.94756982      a Ba 0,50
## Ba 0,75 1.71926315     ab Ba 0,75
## Ba 1    1.62654929    abc    Ba 1
## Ba 10   1.58719230    abc   Ba 10
## Ba 7,5  1.52082616    abc  Ba 7,5
## Ba 2,5  1.49373647    abc  Ba 2,5
## Ba 0,25 1.05502317     bc Ba 0,25
## WC      0.82417544      c      WC
## YC      0.01860719      d      YC

… Et la visualisation graphique :

ggplot(data = cm3, mapping = aes(x = regime, y = indbic)) +
  geom_bar(stat = "identity", color = "blue", fill = "grey", width = 0.6) +
  geom_text(aes(label = groups), vjust = -0.5, size = 4) +
  #ylim(0, 90) +
  xlab("Régimes") + ylab("Lo(Egg shell yellownish index)") +
  theme(axis.text.x = element_text(angle = 45, color = "bbick", vjust = 1, hjust = 1)) +
  theme_bw()

11.5.4 Séance 4

bic4 <- bic %>% filter(seance == "seance 4")

(bic4_out <- bic4 %>% identify_outliers(indbic))

## [1] id         seance     regime     no_oeuf    indbic     is.outlier is.extreme
## <0 rows> (or 0-length row.names)

11.5.4.1 Le modèle

lm4 <- lm(indbic ~ regime, data = bic4)
Anova(lm4)

## Anova Table (Type II tests)
## 
## Response: indbic
##            Sum Sq Df F value Pr(>F)
## regime     1.5086  9  0.2721 0.9752
## Residuals 12.3215 20

La p-value > 0.05 => Pas de différence significative entre les effets des régimes.

shapiro_test(residuals(lm4))

## # A tibble: 1 x 3
##   variable       statistic p.value
##   <chr>              <dbl>   <dbl>
## 1 residuals(lm4)     0.972   0.600

=> Normalité Okay.

11.5.4.2 Comparaisons par paires

cm4 <- (SNK.test(lm4, "regime", group = TRUE))$groups %>% 
  mutate(regime = rownames(.)) %>% 
  select(regime, indbic, groups) %>% 
  as_tibble()
cm4

## # A tibble: 10 x 3
##    regime  indbic groups
##    <chr>    <dbl> <chr> 
##  1 Ba 5     1.56  a     
##  2 WC       1.52  a     
##  3 Ba 0,25  1.45  a     
##  4 Ba 0,50  1.34  a     
##  5 Ba 0,75  1.31  a     
##  6 Ba 1     1.25  a     
##  7 Ba 2,5   1.21  a     
##  8 Ba 10    1.17  a     
##  9 Ba 7,5   0.903 a     
## 10 YC       0.853 a

Visualisation des groupes

ggplot(data = cm4, mapping = aes(x = regime, y = indbic)) +
  geom_bar(stat = "identity", color = "blue", fill = "grey", width = 0.6) +
  #ylim(0, 2) +
  geom_text(aes(label = groups), vjust = -0.5, size = 4) +
  xlab("Régimes") + ylab("Egg shell yellownish index") +
  theme(axis.text.x = element_text(angle = 45, color = "bbick", vjust = 1, hjust = 1)) +
  theme_bw()

11.5.5 Séance 5

bic5 <- bic %>% filter(seance == "seance 5")

(bic5_out <- bic5 %>% identify_outliers(indbic))

## # A tibble: 3 x 7
##   id    seance   regime  no_oeuf indbic is.outlier is.extreme
##   <fct> <fct>    <fct>     <dbl>  <dbl> <lgl>      <lgl>     
## 1 5     seance 5 Ba 0,50       2   5.2  TRUE       FALSE     
## 2 15    seance 5 Ba 2,5        3   4.98 TRUE       FALSE     
## 3 21    seance 5 Ba 7,5        3   5.91 TRUE       FALSE

=> Pas d’observation extrème

#bic5 <- bic5 %>% filter(id != ...)

11.5.5.1 Le modèle

lm5 <- lm(indbic ~ regime, data = bic5)
Anova(lm5)

## Anova Table (Type II tests)
## 
## Response: indbic
##           Sum Sq Df F value Pr(>F)
## regime    31.247  9   1.491  0.218
## Residuals 46.572 20

La p-value est > 0.05 => Pas de différence significative entre les effetes des régimes.

shapiro_test(residuals(lm5))

## # A tibble: 1 x 3
##   variable       statistic p.value
##   <chr>              <dbl>   <dbl>
## 1 residuals(lm5)     0.961   0.323

=> Normalité Okay.

11.5.5.2 Comparaisons par paires

cm5 <- (SNK.test(lm5, "regime", group = TRUE))$groups %>% 
  mutate(regime = rownames(.)) %>% 
  select(regime, indbic, groups) %>% 
  as_tibble()
cm5

## # A tibble: 10 x 3
##    regime  indbic groups
##    <chr>    <dbl> <chr> 
##  1 Ba 0,50   3.35 a     
##  2 Ba 7,5    3.31 a     
##  3 Ba 5      2.74 a     
##  4 Ba 2,5    2.52 a     
##  5 Ba 0,75   1.62 a     
##  6 Ba 0,25   1.26 a     
##  7 Ba 1      1.21 a     
##  8 WC        1.19 a     
##  9 Ba 10     0.62 a     
## 10 YC        0.42 a

11.6 Évolution de l’indice B de l’intérieur de la coquille par régime au cours du temps

11.6.1 Sommaire

bic_ic <- summarySE(bic, 
                   measurevar = "indbic", 
                   groupvars = c("seance", "regime"),
                   na.rm = TRUE)

bic_ic

##      seance  regime N     indbic          sd          se         ci
## 1  seance 1 Ba 0,25 3 13.0100000  1.04971425  0.60605280  2.6076347
## 2  seance 1 Ba 0,50 3 13.8400000  0.94318609  0.54454874  2.3430041
## 3  seance 1 Ba 0,75 3 12.3750000  0.30500000  0.17609183  0.7576620
## 4  seance 1    Ba 1 3 14.1333333  3.22199214  1.86021803  8.0038722
## 5  seance 1   Ba 10 3 13.1000000  0.00000000  0.00000000  0.0000000
## 6  seance 1  Ba 2,5 3 16.0933333  1.08426626  0.62600142  2.6934667
## 7  seance 1    Ba 5 3 15.2300000  0.34770677  0.20074860  0.8637515
## 8  seance 1  Ba 7,5 3 15.2533333  3.18354101  1.83801826  7.9083543
## 9  seance 1      WC 3 40.2100000 26.97910488 15.57639346 67.0198119
## 10 seance 1      YC 3 18.7266667  3.59903228  2.07790225  8.9404918
## 11 seance 2 Ba 0,25 3  1.9366667  0.92088725  0.53167450  2.2876108
## 12 seance 2 Ba 0,50 3  4.2466667  1.32447474  0.76468585  3.2901776
## 13 seance 2 Ba 0,75 3  2.1433333  0.25006666  0.14437605  0.6212000
## 14 seance 2    Ba 1 3  1.8966667  0.20599353  0.11893042  0.5117163
## 15 seance 2   Ba 10 3  2.0600000  0.53018865  0.30610456  1.3170616
## 16 seance 2  Ba 2,5 3  2.1333333  1.12077354  0.64707891  2.7841558
## 17 seance 2    Ba 5 3  3.0133333  0.18175075  0.10493384  0.4514939
## 18 seance 2  Ba 7,5 3  3.5433333  0.12096832  0.06984109  0.3005020
## 19 seance 2      WC 3  3.7566667  3.62108179  2.09063255  8.9952658
## 20 seance 2      YC 3  3.0100000  0.28930952  0.16703293  0.7186847
## 21 seance 3 Ba 0,25 3  3.2950000  1.91500000  1.10562577  4.7571237
## 22 seance 3 Ba 0,50 3  7.3800000  2.87156752  1.65790028  7.1333692
## 23 seance 3 Ba 0,75 3  5.7466667  1.59544142  0.92112853  3.9632962
## 24 seance 3    Ba 1 3  5.5966667  3.11146804  1.79640691  7.7293151
## 25 seance 3   Ba 10 3  4.8900000  0.00000000  0.00000000  0.0000000
## 26 seance 3  Ba 2,5 3  4.4600000  0.29000000  0.16743158  0.7203999
## 27 seance 3    Ba 5 3  7.6000000  1.87000000  1.07964500  4.6453375
## 28 seance 3  Ba 7,5 3  4.6366667  0.88115455  0.50873481  2.1889092
## 29 seance 3      WC 3  2.2800000  0.00000000  0.00000000  0.0000000
## 30 seance 3      YC 3  1.0233333  0.11547005  0.06666667  0.2868435
## 31 seance 4 Ba 0,25 3  1.4466667  0.72154926  0.41658666  1.7924277
## 32 seance 4 Ba 0,50 3  1.3366667  0.20033306  0.11566234  0.4976549
## 33 seance 4 Ba 0,75 3  1.3133333  0.56092186  0.32384839  1.3934071
## 34 seance 4    Ba 1 3  1.2500000  0.93578844  0.54027771  2.3246273
## 35 seance 4   Ba 10 3  1.1666667  0.37859389  0.21858128  0.9404794
## 36 seance 4  Ba 2,5 3  1.2133333  1.26871326  0.73249194  3.1516585
## 37 seance 4    Ba 5 3  1.5600000  0.00000000  0.00000000  0.0000000
## 38 seance 4  Ba 7,5 3  0.9033333  1.02050641  0.58918965  2.5350785
## 39 seance 4      WC 3  1.5166667  1.26894970  0.73262845  3.1522458
## 40 seance 4      YC 3  0.8533333  0.07094599  0.04096069  0.1762396
## 41 seance 5 Ba 0,25 3  1.2633333  1.35433871  0.78192782  3.3643639
## 42 seance 5 Ba 0,50 3  3.3466667  2.10374270  1.21459641  5.2259866
## 43 seance 5 Ba 0,75 3  1.6200000  0.52163205  0.30116441  1.2958059
## 44 seance 5    Ba 1 3  1.2066667  0.59752266  0.34497987  1.4843286
## 45 seance 5   Ba 10 3  0.6200000  0.00000000  0.00000000  0.0000000
## 46 seance 5  Ba 2,5 3  2.5200000  2.17142810  1.25367460  5.3941264
## 47 seance 5    Ba 5 3  2.7400000  1.65737745  0.95688731  4.1171538
## 48 seance 5  Ba 7,5 3  3.3133333  2.91513865  1.68305608  7.2416058
## 49 seance 5      WC 3  1.1933333  0.65317175  0.37710888  1.6225686
## 50 seance 5      YC 3  0.4200000  0.10000000  0.05773503  0.2484138

11.6.2 Visualisation

ggplot(bic_ic, aes(x = seance, y = indbic, colour = regime, group = regime)) + 
  geom_line(size = 1) +
  geom_point(size = 2) +
  ylab("Egg shell yellownish index") +
  theme_bw()

Tendance globale de décroissance d’une séance à l’autre. Vous jugerez.

Nous savons par les analyses pour chaque séance plus haut, que

séance 1 : il existe des différences d’effet entre les régimes
séance 2 : il existe des différences d’effet entre les régimes
séance 3 : il existe des différences d’effet entre les régimes
séance 4 : pas de différences signicatives d’effet entre les régimes
séance 5 : pas de différences signicatives d’effet entre les régimes

Puisque les données ne répondent pas aux conditions pour évaluer les effets des régimes au cours du temps, on négligera l’effet des régimes pour évaluer globalement l’effet du temps sur cet indice.

On pourrait se demander si les indices B mesurés sur l’ensemble des sujets sont significativement différents d’une séance à l’autre (c’est-à-dire avec le temps).

11.6.3 Effet du temps

11.6.3.1 boxplots, facteur temps

bxp <- ggplot(bic, aes(x = seance, y = indbic)) +
  geom_boxplot()
bxp

11.6.3.2 Valeurs aberrantes, facteur temps

bic <- bic %>% mutate(id2 = 1:nrow(.), .before = 1)

bic_out <- bic %>%
  group_by(seance) %>%
  identify_outliers(indbic)
bic_out

## # A tibble: 10 x 8
##    seance     id2 id    regime  no_oeuf indbic is.outlier is.extreme
##    <fct>    <int> <fct> <fct>     <dbl>  <dbl> <lgl>      <lgl>     
##  1 seance 1    29 29    WC            2  30.0  TRUE       TRUE      
##  2 seance 1    30 30    WC            3  70.8  TRUE       TRUE      
##  3 seance 2    35 5     Ba 0,50       2   5.71 TRUE       FALSE     
##  4 seance 2    59 29    WC            2   7.92 TRUE       TRUE      
##  5 seance 3    65 5     Ba 0,50       2  10.4  TRUE       FALSE     
##  6 seance 3    71 11    Ba 1          2   9.15 TRUE       FALSE     
##  7 seance 3    77 17    Ba 5          2   9.47 TRUE       FALSE     
##  8 seance 5   125 5     Ba 0,50       2   5.2  TRUE       FALSE     
##  9 seance 5   135 15    Ba 2,5        3   4.98 TRUE       FALSE     
## 10 seance 5   141 21    Ba 7,5        3   5.91 TRUE       FALSE

=> 2 observations aberrantes extrêmes. Non exclues.

#bic <- bic %>% filter(id2 != c(29, 59))

11.6.3.3 Homogénéité des variances et ANOVA, facteur temps

Les autres conditions ont déjà été vérifiées.

lm <- anova_test(data = bic,
                 dv = indbic,
                 wid = id,
                 within = seance)

get_anova_table(lm)

## ANOVA Table (type III tests)
## 
##   Effect  DFn   DFd      F        p p<.05   ges
## 1 seance 1.13 32.77 50.467 1.51e-08     * 0.588

=> C’est la p-value qui nous intéresse et elle est < 0.05 => Différence significative entre certaines séances.

11.6.3.4 Comparaisons par paires, facteur temps

tph <- bic %>%
  pairwise_t_test(indbic ~ seance, 
                  paired = TRUE,
                  p.adjust.method = "bonferroni")

tph %>% 
  select(group1, group2, p, p.adj, p.adj.signif)

## # A tibble: 10 x 5
##    group1   group2               p        p.adj p.adj.signif
##    <chr>    <chr>            <dbl>        <dbl> <chr>       
##  1 seance 1 seance 2 0.000000054   0.00000054   ****        
##  2 seance 1 seance 3 0.00000278    0.0000278    ****        
##  3 seance 1 seance 4 0.00000000598 0.0000000598 ****        
##  4 seance 1 seance 5 0.0000000225  0.000000225  ****        
##  5 seance 2 seance 3 0.000475      0.005        **          
##  6 seance 2 seance 4 0.0000245     0.000245     ***         
##  7 seance 2 seance 5 0.014         0.138        ns          
##  8 seance 3 seance 4 0.0000000255  0.000000255  ****        
##  9 seance 3 seance 5 0.0000000497  0.000000497  ****        
## 10 seance 4 seance 5 0.087         0.873        ns

11.6.3.5 Boxplots avec p-values

tph <- tph %>% add_xy_position(x = "seance")

ggboxplot(bic, x = "seance", y = "indbic") + 
  #ylim(0, 100) +
  stat_pvalue_manual(tph) +
  labs(subtitle = get_test_label(lm, detailed = TRUE), 
       caption = get_pwc_label(tph))

Chapitre 11 Indice B de l’intérieur de la coquille