Chapitre 5 Masse 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 les masses des coquilles d’oeuf mesurées en 5 séances. Mêmes traitements (régimes).
La question est de savoir si les différents régimes induisent des masses de coquille d’oeuf significativement différentes avec le temps. Mais chaque traitement n’ayant pas été appliqué sur chaque groupe 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 en fonction du temps.
5.1 Les données
mc <- read_csv("data/oeuf.csv")
mc <- mc %>%
select(seance, regime, no_oeuf, masse_coq) %>%
mutate(id = rep(1:30, 5), .before = 1,
id = factor(id),
seance = factor(seance),
regime = factor(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.
## Rows: 150
## Columns: 5
## $ id <fct> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1~
## $ seance <fct> seance 1, seance 1, seance 1, seance 1, seance 1, seance 1, ~
## $ regime <fct> "Ba 0,25", "Ba 0,25", "Ba 0,25", "Ba 0,50", "Ba 0,50", "Ba 0~
## $ no_oeuf <dbl> 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, ~
## $ masse_coq <dbl> 1.187, 1.313, 1.431, 1.317, 1.176, 1.463, 1.107, 1.054, 1.08~5.2 Visualisation boxplots
bxp <- ggplot(mc, aes(x = regime, y = masse_coq, fill = regime)) +
geom_boxplot() +
facet_grid(seance ~ .) +
theme(axis.text.x = element_text(angle = 90, color = "black", vjust = 0.5, hjust = 1)) +
theme_bw()
bxp
Il y a des variations notables d’une séance à l’autre pour certains régimes.
5.3 Détection des valeurs aberrantes extrêmes
## [1] seance regime id no_oeuf masse_coq is.outlier is.extreme
## <0 rows> (or 0-length row.names)=> Pas de valeurs aberrantes extrêmes pour toutes les séances.
5.4 Conditions de l’ANOVA
5.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.
## # A tibble: 5 x 4
## seance variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 seance 1 masse_coq 0.962 0.356
## 2 seance 2 masse_coq 0.964 0.385
## 3 seance 3 masse_coq 0.954 0.214
## 4 seance 4 masse_coq 0.914 0.0193
## 5 seance 5 masse_coq 0.972 0.606=> Normalité confirmée pour toutes les séances sauf la 4.
Créer des QQ-plots pour chaque point par séance
Tous les points se situent approximativement le long de la ligne de référence. Mais on regardera de plus près la séance 4 au moment venu.
5.4.2 Homogénéité des variances
## # A tibble: 5 x 5
## seance df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 seance 1 9 20 1.16 0.373
## 2 seance 2 9 20 0.408 0.916
## 3 seance 3 9 20 0.988 0.480
## 4 seance 4 9 20 0.565 0.810
## 5 seance 5 9 20 0.658 0.736=> Toutes les valeurs p sont > 0.05 => toutes les variances sont homogènes.
Les conditions de la validité d’une ANOVA étant remplies, les interprétations seront donc valides. Le cas de la séance 4 vera une vérification supplémentaire.
5.5 ANOVA à 1 facteur séance par séance
5.5.1 Séance 1
5.5.1.1 Le modèle
## Anova Table (Type II tests)
##
## Response: masse_coq
## Sum Sq Df F value Pr(>F)
## regime 0.47716 9 5.306 0.0009252 ***
## Residuals 0.19984 20
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La p-value < 0.01 => différence très significative entre les effets d’au moins 2 régimes sur la masse de la coquille à la séance 1.
5.5.1.2 Comparaisons par paires
Comparaisons des moyennes par paires, Student - Newman - Keuls.
cm1 <- (SNK.test(lm1, "regime", group = TRUE))$groups %>%
mutate(regime = rownames(.)) %>%
select(regime, masse_coq, groups) %>%
as_tibble()
cm1## # A tibble: 10 x 3
## regime masse_coq groups
## <chr> <dbl> <chr>
## 1 YC 1.49 a
## 2 WC 1.42 a
## 3 Ba 10 1.38 ab
## 4 Ba 0,50 1.32 abc
## 5 Ba 0,25 1.31 abc
## 6 Ba 1 1.26 abc
## 7 Ba 5 1.25 abc
## 8 Ba 2,5 1.15 bc
## 9 Ba 7,5 1.12 bc
## 10 Ba 0,75 1.08 c5.5.1.3 Visualisation des groupes, bareplots avec labels
Figure pas nécessaire.
ggplot(data = cm1, mapping = aes(x = regime, y = masse_coq)) +
geom_bar(stat = "identity", color = "blue", fill = "grey", width = 0.6) +
geom_text(aes(label = groups), vjust = -0.5, size = 4) +
ylim(0, 1.6) +
xlab("Régimes") + ylab("Masse de la coquille (g)") +
theme(axis.text.x = element_text(angle = 45, color = "black", vjust = 1, hjust = 1)) +
theme_bw()
5.5.2 Séance 2
5.5.2.1 Le modèle
## Anova Table (Type II tests)
##
## Response: masse_coq
## Sum Sq Df F value Pr(>F)
## regime 0.35177 9 2.582 0.03707 *
## Residuals 0.30275 20
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La p-value est < 0.05 => Différence significative entre les effetes d’au moins 2 régimes.
5.5.2.2 Comparaisons par paires
cm2 <- (SNK.test(lm2, "regime", group = TRUE))$groups %>%
mutate(regime = rownames(.)) %>%
select(regime, masse_coq, groups) %>%
as_tibble()
cm2## # A tibble: 10 x 3
## regime masse_coq groups
## <chr> <dbl> <chr>
## 1 YC 1.56 a
## 2 Ba 1 1.48 ab
## 3 Ba 0,25 1.41 ab
## 4 Ba 5 1.38 ab
## 5 Ba 0,75 1.37 ab
## 6 WC 1.30 ab
## 7 Ba 7,5 1.3 ab
## 8 Ba 10 1.28 ab
## 9 Ba 0,50 1.22 ab
## 10 Ba 2,5 1.20 b5.5.2.3 Visualisation des groupes
ggplot(data = cm2, mapping = aes(x = regime, y = masse_coq)) +
geom_bar(stat = "identity", color = "blue", fill = "grey", width = 0.6) +
ylim(0, 1.6) +
geom_text(aes(label = groups), vjust = -0.5, size = 4) +
xlab("Régimes") + ylab("Masse de la coquille (g)") +
theme(axis.text.x = element_text(angle = 45, color = "black", vjust = 1, hjust = 1)) +
theme_bw()
5.5.3 Séance 3
5.5.3.1 Le modèle
## Anova Table (Type II tests)
##
## Response: masse_coq
## Sum Sq Df F value Pr(>F)
## regime 0.35845 9 3.1607 0.01537 *
## Residuals 0.25202 20
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La p-value < 0.05 => Différence significative entre les effetes d’au moins 2 régimes.
5.5.3.2 Comparaisons par paires
cm3 <- (SNK.test(lm3, "regime", group = TRUE))$groups %>%
mutate(regime = rownames(.)) %>%
select(regime, masse_coq, groups) %>%
as_tibble()
cm3## # A tibble: 10 x 3
## regime masse_coq groups
## <chr> <dbl> <chr>
## 1 Ba 2,5 1.53 a
## 2 Ba 0,50 1.42 ab
## 3 Ba 1 1.39 ab
## 4 Ba 0,25 1.39 ab
## 5 Ba 7,5 1.37 ab
## 6 WC 1.33 ab
## 7 Ba 5 1.31 ab
## 8 YC 1.24 ab
## 9 Ba 0,75 1.22 ab
## 10 Ba 10 1.13 b5.5.3.3 Visualisation des groupes
ggplot(data = cm3, mapping = aes(x = regime, y = masse_coq)) +
geom_bar(stat = "identity", color = "blue", fill = "grey", width = 0.6) +
geom_text(aes(label = groups), vjust = -0.5, size = 4) +
ylim(0, 1.6) +
xlab("Régimes") + ylab("Masse de la coquille (g)") +
theme(axis.text.x = element_text(angle = 45, color = "black", vjust = 1, hjust = 1)) +
theme_bw()
5.5.4 Séance 4
## # A tibble: 1 x 7
## id seance regime no_oeuf masse_coq is.outlier is.extreme
## <fct> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 3 seance 4 Ba 0,25 3 1.65 TRUE FALSEExcluons cette observation pour apprécier.
## # A tibble: 1 x 3
## variable statistic p
## <chr> <dbl> <dbl>
## 1 masse_coq 0.940 0.0981=> Normalité OK !
## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 9 19 0.719 0.686=> Homogénéité des variances toujours OK !
5.5.4.1 Le modèle
## Anova Table (Type II tests)
##
## Response: masse_coq
## Sum Sq Df F value Pr(>F)
## regime 0.23631 9 1.2033 0.3485
## Residuals 0.41461 19La p-value > 0.05 => différence pas significative entre les effets des régimes.
5.5.4.2 Comparaisons par paires, séance 4
cm4 <- (SNK.test(lm4, "regime", group = TRUE))$groups %>%
mutate(regime = rownames(.)) %>%
select(regime, masse_coq, groups) %>%
as_tibble()
cm4## # A tibble: 10 x 3
## regime masse_coq groups
## <chr> <dbl> <chr>
## 1 Ba 1 1.30 a
## 2 Ba 0,25 1.30 a
## 3 Ba 2,5 1.27 a
## 4 Ba 7,5 1.24 a
## 5 Ba 10 1.22 a
## 6 YC 1.18 a
## 7 WC 1.13 a
## 8 Ba 0,50 1.08 a
## 9 Ba 0,75 1.07 a
## 10 Ba 5 1.04 a5.5.4.3 Visualisation des groupes
ggplot(data = cm4, mapping = aes(x = regime, y = masse_coq)) +
geom_bar(stat = "identity", color = "blue", fill = "grey", width = 0.6) +
ylim(0, 1.5) +
geom_text(aes(label = groups), vjust = -0.5, size = 4) +
xlab("Régimes") + ylab("Masse de la coquille (g)") +
theme(axis.text.x = element_text(angle = 45, color = "black", vjust = 1, hjust = 1)) +
theme_bw()
5.5.5 Séance 5
5.5.5.1 Le modèle
## Anova Table (Type II tests)
##
## Response: masse_coq
## Sum Sq Df F value Pr(>F)
## regime 0.14107 9 1.1202 0.3935
## Residuals 0.27985 20La p-value > 0.05 => différence pas significative entre les effets des régimes.
5.5.5.2 Comparaisons par paires
cm5 <- (SNK.test(lm5, "regime", group = TRUE))$groups %>%
mutate(regime = rownames(.)) %>%
select(regime, masse_coq, groups) %>%
as_tibble()
cm5## # A tibble: 10 x 3
## regime masse_coq groups
## <chr> <dbl> <chr>
## 1 Ba 0,50 1.51 a
## 2 Ba 10 1.48 a
## 3 Ba 1 1.42 a
## 4 WC 1.40 a
## 5 Ba 0,75 1.39 a
## 6 Ba 7,5 1.37 a
## 7 YC 1.33 a
## 8 Ba 0,25 1.33 a
## 9 Ba 2,5 1.30 a
## 10 Ba 5 1.30 a5.5.5.3 Visualisation
ggplot(data = cm5, mapping = aes(x = regime, y = masse_coq)) +
geom_bar(stat = "identity", color = "blue", fill = "grey", width = 0.6) +
ylim(0, 1.6) +
geom_text(aes(label = groups), vjust = -0.5, size = 4) +
xlab("Régimes") + ylab("Masse de la coquille (g)") +
theme(axis.text.x = element_text(angle = 45, color = "black", vjust = 1, hjust = 1)) +
theme_bw()
5.6 Évolution de la masse de la coquille par régime au cours du temps
5.6.1 Sommaire
mc_ic <- summarySE(mc,
measurevar = "masse_coq",
groupvars = c("seance", "regime"),
na.rm = TRUE)
mc_ic## seance regime N masse_coq sd se ci
## 1 seance 1 Ba 0,25 3 1.310333 0.122021856 0.0704493514 0.303119094
## 2 seance 1 Ba 0,50 3 1.318667 0.143507259 0.0828539545 0.356491794
## 3 seance 1 Ba 0,75 3 1.080667 0.026501572 0.0153006899 0.065833555
## 4 seance 1 Ba 1 3 1.263000 0.128175661 0.0740022522 0.318405993
## 5 seance 1 Ba 10 3 1.381000 0.000000000 0.0000000000 0.000000000
## 6 seance 1 Ba 2,5 3 1.148333 0.106326541 0.0613876571 0.264129770
## 7 seance 1 Ba 5 3 1.255000 0.037643060 0.0217332311 0.093510546
## 8 seance 1 Ba 7,5 3 1.121000 0.026851443 0.0155026879 0.066702683
## 9 seance 1 WC 3 1.420667 0.110028784 0.0635251481 0.273326652
## 10 seance 1 YC 3 1.488333 0.147500282 0.0851593278 0.366411014
## 11 seance 2 Ba 0,25 3 1.406333 0.063877487 0.0368796843 0.158680474
## 12 seance 2 Ba 0,50 3 1.218000 0.099954990 0.0577090403 0.248301960
## 13 seance 2 Ba 0,75 3 1.368667 0.095311769 0.0550282756 0.236767560
## 14 seance 2 Ba 1 3 1.480667 0.134838175 0.0778488564 0.334956595
## 15 seance 2 Ba 10 3 1.281000 0.109672239 0.0633192967 0.272440945
## 16 seance 2 Ba 2,5 3 1.195333 0.065033325 0.0375470076 0.161551735
## 17 seance 2 Ba 5 3 1.382000 0.138683092 0.0800687205 0.344507899
## 18 seance 2 Ba 7,5 3 1.300000 0.163450298 0.0943680737 0.406033050
## 19 seance 2 WC 3 1.304333 0.204756766 0.1182163741 0.508644005
## 20 seance 2 YC 3 1.562000 0.076863515 0.0443771713 0.190939557
## 21 seance 3 Ba 0,25 3 1.389000 0.063174362 0.0364737348 0.156933815
## 22 seance 3 Ba 0,50 3 1.415333 0.134440817 0.0776194420 0.333969504
## 23 seance 3 Ba 0,75 3 1.219000 0.057087652 0.0329595712 0.141813589
## 24 seance 3 Ba 1 3 1.394667 0.135692053 0.0783418435 0.337077747
## 25 seance 3 Ba 10 3 1.127000 0.000000000 0.0000000000 0.000000000
## 26 seance 3 Ba 2,5 3 1.534000 0.121000000 0.0698593826 0.300580663
## 27 seance 3 Ba 5 3 1.314500 0.000500000 0.0002886751 0.001242069
## 28 seance 3 Ba 7,5 3 1.374667 0.240803516 0.1390279748 0.598189095
## 29 seance 3 WC 3 1.330000 0.000000000 0.0000000000 0.000000000
## 30 seance 3 YC 3 1.244333 0.098215749 0.0567048891 0.243981446
## 31 seance 4 Ba 0,25 3 1.413333 0.202376711 0.1168422488 0.502731621
## 32 seance 4 Ba 0,50 3 1.078000 0.160411970 0.0926138938 0.398485423
## 33 seance 4 Ba 0,75 3 1.068333 0.121911170 0.0703854467 0.302844134
## 34 seance 4 Ba 1 3 1.297333 0.055365453 0.0319652589 0.137535409
## 35 seance 4 Ba 10 3 1.216667 0.125707332 0.0725771620 0.312274324
## 36 seance 4 Ba 2,5 3 1.266667 0.141923688 0.0819396797 0.352557987
## 37 seance 4 Ba 5 3 1.045000 0.000000000 0.0000000000 0.000000000
## 38 seance 4 Ba 7,5 3 1.241333 0.342077964 0.1974988045 0.849768770
## 39 seance 4 WC 3 1.132333 0.008082904 0.0046666667 0.020079046
## 40 seance 4 YC 3 1.177667 0.103001618 0.0594680119 0.255870204
## 41 seance 5 Ba 0,25 3 1.326667 0.095552778 0.0551674219 0.237366259
## 42 seance 5 Ba 0,50 3 1.509667 0.109235220 0.0630669838 0.271355330
## 43 seance 5 Ba 0,75 3 1.391000 0.056204982 0.0324499615 0.139620915
## 44 seance 5 Ba 1 3 1.415667 0.069176104 0.0399388421 0.171842968
## 45 seance 5 Ba 10 3 1.480000 0.000000000 0.0000000000 0.000000000
## 46 seance 5 Ba 2,5 3 1.302333 0.113058982 0.0652746335 0.280854080
## 47 seance 5 Ba 5 3 1.299333 0.138911243 0.0802004433 0.345074656
## 48 seance 5 Ba 7,5 3 1.370667 0.273234942 0.1577522671 0.678753223
## 49 seance 5 WC 3 1.404667 0.048273526 0.0278707334 0.119918087
## 50 seance 5 YC 3 1.327000 0.043000000 0.0248260616 0.1068179225.6.2 Visualisation
ggplot(mc_ic, aes(x = seance, y = masse_coq, colour = regime, group = regime)) +
geom_line(size = 1) +
geom_point(size = 2) +
ylab("Masse de la coquille (g)") +
theme_bw()
Il ne semble pas se dégager une tendance concrète entre les différentes séances, ou difficile à exploiter. 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 (1 valeur aberrante exclue)
- séance 5 : pas de différences signicatives
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 les masses des coquilles d’oeuf.
On pourrait se demander si les masses des coquilles mesurées sur l’ensemble des sujets sont significativement différentes d’une séance à l’autre (c’est-à-dire avec le temps).
5.6.3 Effet du temps
5.6.3.1 boxplots, facteur temps
5.6.3.2 Valeurs aberrantes, facteur temps
mc <- mc %>% mutate(id2 = 1:nrow(.), .before = 1)
mc_out <- mc %>%
group_by(seance) %>%
identify_outliers(masse_coq)
mc_out## # A tibble: 3 x 8
## seance id2 id regime no_oeuf masse_coq is.outlier is.extreme
## <fct> <int> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 seance 3 73 13 Ba 2,5 1 1.66 TRUE FALSE
## 2 seance 4 93 3 Ba 0,25 3 1.65 TRUE FALSE
## 3 seance 5 140 20 Ba 7,5 2 1.68 TRUE FALSE=> 3 observations aberrantes mais pas extrêmes pour les séances 3, 4 et 5 ; une par séance. On pourrait les exclure. Mais je les ai conservé pour la suite.
5.6.3.3 Homogénéité des variances et ANOVA, facteur temps
Les autres conditions ont déjà été vérifiées. La fonction anova_test() réalise également le test de sphéricité de Mauchly.
lm <- anova_test(data = mc,
dv = masse_coq, # dependant variable, num
wid = id, # identificateur de cas/échantillon (facteur)
within = seance) # facteur de groupement intra-sujets
get_anova_table(lm)## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 seance 4 116 7.6 1.79e-05 * 0.169=> C’est la p-value qui nous intéresse et elle est < 0.05 => différence significative entre certaines séances.
5.6.3.4 Comparaisons par paires, facteur temps
tph <- mc %>%
pairwise_t_test(masse_coq ~ 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.066 0.663 ns
## 2 seance 1 seance 3 0.196 1 ns
## 3 seance 1 seance 4 0.042 0.422 ns
## 4 seance 1 seance 5 0.004 0.041 *
## 5 seance 2 seance 3 0.698 1 ns
## 6 seance 2 seance 4 0.000175 0.002 **
## 7 seance 2 seance 5 0.404 1 ns
## 8 seance 3 seance 4 0.000297 0.003 **
## 9 seance 3 seance 5 0.173 1 ns
## 10 seance 4 seance 5 0.00011 0.001 **5.6.3.5 Boxplots avec p-values
tph <- tph %>% add_xy_position(x = "seance")
ggboxplot(mc, x = "seance", y = "masse_coq") +
stat_pvalue_manual(tph) +
labs(subtitle = get_test_label(lm, detailed = TRUE),
caption = get_pwc_label(tph))
Un peu difficile à exploiter mais on pourrait se baser sur le tableau du code ci-dessus pour les déméler …