Document complémentaire au module 9 du cours SDD I de 2025-2026. Distribué sous licence CC BY-NC-SA 4.0.

Veuillez vous référer au cours en ligne pour les explications et les interprétations de cette analyse.

Installer un environnement R adéquat pour reproduire cette analyse.

ANOVA à un facteur

# Configuration de R pour le dialecte SciViews::R avec les modules infer et model
SciViews::R("infer", "model", lang = "fr")

# Lecture du jeu de données depuis le package MASS
crabs <- read("crabs", package = "MASS", lang = "fr")
# Création d'une variable de groupe combinant espèce et sexe avec label correct
crabs2 <- smutate(crabs, group  = labelise(
    factor(paste(species, sex, sep = "-")),
    "Groupe espèce - sexe", units = NA))

# Graphique de comparaison de la largeur à l'arrière en fonction du groupe
chart(data = crabs2, rear ~ group) +
  geom_violin() +
  geom_jitter(width = 0.05, alpha = 0.5) +
  stat_summary(geom = "point", fun = "mean", color = "red", size = 3)

# Calcul d'une variable aspect (ratio largeur arrière / largeur max) avec label
crabs2 %>.%
  smutate(., aspect = labelise(
    as.numeric(rear / width),
    "Ratio largeur arrière / max", units = NA)) %>.%
  sselect(., species, sex, group, aspect) ->
  crabs2
# Exploration des données
skimr::skim(crabs2)
Data summary
Name crabs2
Number of rows 200
Number of columns 4
_______________________
Column type frequency:
factor 3
numeric 1
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
species 0 1 FALSE 2 B: 100, O: 100
sex 0 1 FALSE 2 F: 100, M: 100
group 0 1 FALSE 4 B-F: 50, B-M: 50, O-F: 50, O-M: 50

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
aspect 0 1 0.35 0.03 0.28 0.32 0.36 0.38 0.41 ▂▅▃▇▂ 
# Graphique du ratio largeur arrière / largeur max en fonction du groupe
chart(data = crabs2, aspect ~ group) +
  geom_violin() +
  geom_jitter(width = 0.05, alpha = 0.5) +
  stat_summary(geom = "point", fun = "mean", color = "red", size = 3)

# Test de Bartlett pour vérifier l'homoscédasticité
bartlett.test(data = crabs2, aspect ~ group)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  aspect by group
## Bartlett's K-squared = 24.532, df = 3, p-value = 1.935e-05
# Nouveau test d'homoscédasticité sur données log-transformées
crabs2 <- smutate(crabs2, log_aspect = log(aspect))
bartlett.test(data = crabs2, log_aspect ~ group)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  log_aspect by group
## Bartlett's K-squared = 37.891, df = 3, p-value = 2.981e-08
# Essai d'une transformation puissance 5 pour stabiliser la variance
crabs2 <- smutate(crabs2, aspect5 = aspect^5)
bartlett.test(data = crabs2, aspect5 ~ group)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  aspect5 by group
## Bartlett's K-squared = 1.7948, df = 3, p-value = 0.6161
# Graphique des données avec transformation puissance 5 de aspect
chart(data = crabs2, aspect5 ~ group) +
  geom_violin() +
  geom_jitter(width = 0.05, alpha = 0.5) +
  stat_summary(geom = "point", fun = "mean", color = "red", size = 3) +
  ylab("(ratio largeur arrière/max)^5")

# Tableau descriptif des données en vue d'une ANOVA
crabs2 %>.%
  sgroup_by(., group) %>.%
  ssummarise(.,
    mean  = fmean(aspect5),
    sd    = fsd(aspect5),
    count = fsum(!is.na(aspect5)))
## # A data.trame: [4 × 4]
##   group    mean      sd count
##   <fct>   <dbl>   <dbl> <int>
## 1 B-F   0.00727 0.00115    50
## 2 B-M   0.00363 0.00124    50
## 3 O-F   0.00811 0.00138    50
## 4 O-M   0.00427 0.00130    50
# Analyse de variance à un facteur (présentation avec tabularise)
anova(crabs2_anova <- lm(data = crabs2, aspect5 ~ group)) |>
  tabularise()

Analyse de la variance
Réponse : aspect5

Terme

Ddl

Somme des carrés

Carrés moyens

Valeur de Fobs.

Valeur de p

group

3

0.000727

2.42·10-04

151

< 2·10-16

***

Résidus

196

0.000316

1.61·10-06

0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

# Graphique quantile-quantile de distribution des résidus
chart$qqplot(crabs2_anova, lang = "fr")

Tests post-hoc

# Comparaisons multiples avec la méthode HSD de Tukey
# Version textuelle avec summary()
summary(crabs2_posthoc <- confint(multcomp::glht(crabs2_anova,
  linfct = multcomp::mcp(group = "Tukey"))))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = aspect5 ~ group, data = crabs2)
## 
## Linear Hypotheses:
##                  Estimate Std. Error t value Pr(>|t|)    
## B-M - B-F == 0 -0.0036378  0.0002538 -14.331   <0.001 ***
## O-F - B-F == 0  0.0008441  0.0002538   3.326   0.0057 ** 
## O-M - B-F == 0 -0.0029979  0.0002538 -11.810   <0.001 ***
## O-F - B-M == 0  0.0044820  0.0002538  17.657   <0.001 ***
## O-M - B-M == 0  0.0006399  0.0002538   2.521   0.0596 .  
## O-M - O-F == 0 -0.0038420  0.0002538 -15.136   <0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
# Version graphique avec plot()
oma <- par(oma = c(0, 5.1, 0, 0))
plot(crabs2_posthoc)

par(oma)
rm(oma)

Test de Kruskal-Wallis

# Calcul des rangs manuellement (pour illustration)
# Un échantillon exemple au hasard
x <- c(4.5, 2.1, 0.5, 2.4, 2.1, 3.5)
# Tri par ordre croissant
sort(x)
## [1] 0.5 2.1 2.1 2.4 3.5 4.5
# Remplacement par les rangs
rank(sort(x))
## [1] 1.0 2.5 2.5 4.0 5.0 6.0
kruskal.test(data = crabs2, aspect ~ group)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  aspect by group
## Kruskal-Wallis chi-squared = 139.68, df = 3, p-value < 2.2e-16
# Comparaisons multiples non paramétriques après un test de Kruskal-Wallis
summary(crabs2_kw_comp <- nparcomp::nparcomp(data = crabs2, aspect ~ group))
## 
##  #------Nonparametric Multiple Comparisons for relative contrast effects-----# 
##  
##  - Alternative Hypothesis:  True relative contrast effect p is not equal to 1/2 
##  - Type of Contrast : Tukey 
##  - Confidence level: 95 % 
##  - Method = Logit - Transformation 
##  - Estimation Method: Pairwise rankings 
##  
##  #---------------------------Interpretation----------------------------------# 
##  p(a,b) > 1/2 : b tends to be larger than a 
##  #---------------------------------------------------------------------------# 
##  
## 
##  #------------Nonparametric Multiple Comparisons for relative contrast effects----------# 
##  
##  - Alternative Hypothesis:  True relative contrast effect p is not equal to 1/2 
##  - Estimation Method: Global Pseudo ranks 
##  - Type of Contrast : Tukey 
##  - Confidence Level: 95 % 
##  - Method = Logit - Transformation 
##  
##  - Estimation Method: Pairwise rankings 
##  
##  #---------------------------Interpretation--------------------------------------------# 
##  p(a,b) > 1/2 : b tends to be larger than a 
##  #-------------------------------------------------------------------------------------# 
##  
##  #----Data Info-------------------------------------------------------------------------# 
##   Sample Size
## 1    B-F   50
## 2    B-M   50
## 3    O-F   50
## 4    O-M   50
## 
##  #----Contrast--------------------------------------------------------------------------# 
##           B-F B-M O-F O-M
## B-M - B-F  -1   1   0   0
## O-F - B-F  -1   0   1   0
## O-M - B-F  -1   0   0   1
## O-F - B-M   0  -1   1   0
## O-M - B-M   0  -1   0   1
## O-M - O-F   0   0  -1   1
## 
##  #----Analysis--------------------------------------------------------------------------# 
##       Comparison Estimator Lower Upper Statistic      p.Value
## 1 p( B-F , B-M )     0.017 0.004 0.069 -7.229567 1.338263e-12
## 2 p( B-F , O-F )     0.695 0.542 0.815  3.257267 6.331133e-03
## 3 p( B-F , O-M )     0.051 0.018 0.137 -7.014700 7.743917e-12
## 4 p( B-M , O-F )     0.990 0.949 0.998  7.210499 2.363332e-12
## 5 p( B-M , O-M )     0.654 0.501 0.781  2.606445 4.797910e-02
## 6 p( O-F , O-M )     0.027 0.008 0.088 -7.385440 4.904965e-13
## 
##  #----Overall---------------------------------------------------------------------------# 
##   Quantile      p.Value
## 1 2.588746 4.904965e-13
## 
##  #--------------------------------------------------------------------------------------#
plot(crabs2_kw_comp)