SDD I module 10 : Variance II

Document complémentaire au module 10 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 à deux facteurs sans interactions

# Chargement du dialecte SciViews::R avec les modules infer et model
SciViews::R("infer", "model", lang = "fr")
# Lecture et transformation des données (calcul de aspect^5)
read("crabs", package = "MASS", lang = "fr") %>.%
  smutate(., aspect = labelise(
    as.numeric(rear / width),
    "Ratio largeur arrière / max", units = NA)) %>.%
  smutate(., aspect5 = labelise(
    aspect^5,
    "(Ratio largeur arrière /max)^5", units = NA)) %>.%
  sselect(., species, sex, aspect, aspect5) ->
  crabs2

# Graphique de base pour visualiser les interactions
#chart$base(interaction.plot(crabs2$species, crabs2$sex, crabs2$aspect5))

# Version avec chart
crabs2 %>.%
  sgroup_by(., species, sex) %>.%
  ssummarise(., aspect5_groups = mean(aspect5)) %>.%
  print(.) %>.% # Tableau des moyennes par groupes
  chart(data = ., aspect5_groups ~ species %col=% sex %group=% sex) +
    geom_line() +
    geom_point()

## # A data.trame: [4 × 3]
##   species sex   aspect5_groups
##   <fct>   <fct>          <dbl>
## 1 B       F            0.00727
## 2 B       M            0.00363
## 3 O       F            0.00811
## 4 O       M            0.00427

# Résumé des données en vue de réaliser une ANOVA à deux facteurs
crabs2 %>.%
  sgroup_by(., species, sex) %>.%
  ssummarise(.,
    "Moyenne aspect^5"  = fmean(aspect5),
    "Variance aspect^5" = fvar(aspect5),
    "N"                 = fnobs(aspect5))

## # A data.trame: [4 × 5]
##   species sex   `Moyenne aspect^5` `Variance aspect^5`     N
##   <fct>   <fct>              <dbl>               <dbl> <int>
## 1 B       F                0.00727          0.00000132    50
## 2 B       M                0.00363          0.00000153    50
## 3 O       F                0.00811          0.00000192    50
## 4 O       M                0.00427          0.00000168    50

# Visualisation des données en vue d'une ANOVA à deux facteurs
chart(data = crabs2, aspect5 ~ species | sex) +
  geom_boxplot()

# Version améliorée avec les observations et les moyennes
chart(data = crabs2, aspect5 ~ species | sex) +
  geom_boxplot() +
  geom_jitter(width = 0.05, alpha = 0.3) +
  stat_summary(geom = "point", fun = "mean", color = "red", size = 3)

# Test d'homoscédasticité de Batlett pour l'ANOVA à deux facteurs
bartlett.test(data = crabs2, aspect5 ~ interaction(species, sex))

## 
##  Bartlett test of homogeneity of variances
## 
## data:  aspect5 by interaction(species, sex)
## Bartlett's K-squared = 1.7948, df = 3, p-value = 0.6161

# ANOVA à deux facteurs sans interaction
anova(crabs2_anova2 <- lm(data = crabs2, aspect5 ~ species + sex)) %>.%
  tabularise(.)

Analyse de la variance Réponse : aspect5
Terme	Ddl	Somme des carrés	Carrés moyens	Valeur de Fobs.	Valeur de p
species	1	2.75·10-05	2.75·10-05	17.2	5.12·10-05	***
sex	1	6.99·10-04	6.99·10-04	435.7	< 2·10-16	***
Résidus	197	3.16·10-04	1.61·10-06
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05

# Tests de comparaisons multiples selon la méthode HSD de Tukey
summary(crabs2_posthoc2 <- confint(multcomp::glht(crabs2_anova2,
  linfct = multcomp::mcp(species = "Tukey", sex = "Tukey"))))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = aspect5 ~ species + sex, data = crabs2)
## 
## Linear Hypotheses:
##                       Estimate Std. Error t value Pr(>|t|)    
## species: O - B == 0  0.0007420  0.0001792   4.141 0.000102 ***
## sex: M - F == 0     -0.0037399  0.0001792 -20.872  < 1e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

oma <- par(oma = c(0, 5.1, 0, 0))
plot(crabs2_posthoc2)

par(oma)
rm(oma)

# Graphique quantile-quantile de notre ANOVA à deux facteurs sans interactions
chart$qqplot(crabs2_anova2)

ANOVA à deux facteurs croisés complet

# Dénombrement des observations par niveaux (ici le plan est balancé)
scount(crabs2, species, sex)

## # A data.trame: [4 × 3]
##   species sex       n
##   <fct>   <fct> <int>
## 1 B       F        50
## 2 B       M        50
## 3 O       F        50
## 4 O       M        50

# ANOVA à deux facteurs, modèle complet avec interactions
crabs2_anova2comp <- lm(data = crabs2, aspect5 ~ species * sex)
anova(crabs2_anova2comp) %>.%
  tabularise(.)

Analyse de la variance Réponse : aspect5
Terme	Ddl	Somme des carrés	Carrés moyens	Valeur de Fobs.	Valeur de p
species	1	2.75·10-05	2.75·10-05	17.091	5.28·10-05	***
sex	1	6.99·10-04	6.99·10-04	434.161	< 2·10-16	***
species:sex	1	5.21·10-07	5.21·10-07	0.324	5.70·10-01
Résidus	196	3.16·10-04	1.61·10-06
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05

ANOVA à deux facteurs hiérarchisés

# Lecture des données eggs depuis le package faraway
eggs <- read("eggs", package = "faraway")
# Exploration des données
skimr::skim(eggs)

Data summary
Name	eggs
Number of rows	48
Number of columns	4
_______________________
Column type frequency:
factor	3
numeric	1
________________________
Group variables	None

Variable type: factor

skim_variable	complete_rate	ordered	n_unique	top_counts
Lab	1	FALSE	6	I: 8, II: 8, III: 8, IV: 8
Technician	1	FALSE	2	one: 24, two: 24
Sample	1	FALSE	2	G: 24, H: 24

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Fat	0	1	0.39	0.15	0.06	0.31	0.37	0.43	0.8	▂▅▇▂▁

# Correction de l'encodage (techniciens imbriqués dans les laboratoires)
eggs <- smutate(eggs, Technician = interaction(Lab, Technician))
# Vérification des données
skimr::skim(eggs)

Data summary
Name	eggs
Number of rows	48
Number of columns	4
_______________________
Column type frequency:
factor	3
numeric	1
________________________
Group variables	None

Variable type: factor

skim_variable	complete_rate	ordered	n_unique	top_counts
Lab	1	FALSE	6	I: 8, II: 8, III: 8, IV: 8
Technician	1	FALSE	12	I.o: 4, II.: 4, III: 4, IV.: 4
Sample	1	FALSE	2	G: 24, H: 24

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Fat	0	1	0.39	0.15	0.06	0.31	0.37	0.43	0.8	▂▅▇▂▁

# Graphique de description des données adéquat pour une ANOVA à facteurs hiérarchisés
chart(data = eggs, Fat ~ Lab %col=% Technician) +
  geom_jitter(width = 0.05, alpha = 0.5) +
  stat_summary(geom = "point", fun = "mean", color = "red", size = 3)

# Test d'homoscédasticité de Bartlett
bartlett.test(data = eggs, Fat ~ Technician)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Fat by Technician
## Bartlett's K-squared = 13.891, df = 11, p-value = 0.2391

# ANOVA avec facteur technicien imbriqué dans le facteur laboratoire
eggs_anova <- lm(data = eggs, Fat ~ Lab + Technician %in% Lab)
anova(eggs_anova) %>.%
  tabularise(.)

Analyse de la variance Réponse : Fat
Terme	Ddl	Somme des carrés	Carrés moyens	Valeur de Fobs.	Valeur de p
Lab	5	0.443	0.08861	9.59	6.99·10-06	***
Lab:Technician	6	0.247	0.04125	4.46	1.79·10-03	**
Résidus	36	0.333	0.00924
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05

# Graphique quantile-quantile des résidus de notre ANOVA à facteurs hiérarchisés
chart$qqplot(eggs_anova)

# Comparaisons multiples par HSD de Tukey en partant de aov()
eggs_aov <- aov(data = eggs, Fat ~ Lab + Technician %in% Lab)
(eggs_posthoc <- TukeyHSD(eggs_aov, "Lab"))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Fat ~ Lab + Technician %in% Lab, data = eggs)
## 
## $Lab
##            diff         lwr          upr     p adj
## II-I   -0.24000 -0.38459081 -0.095409195 0.0002088
## III-I  -0.17250 -0.31709081 -0.027909195 0.0116356
## IV-I   -0.20375 -0.34834081 -0.059159195 0.0019225
## V-I    -0.22625 -0.37084081 -0.081659195 0.0004902
## VI-I   -0.31250 -0.45709081 -0.167909195 0.0000021
## III-II  0.06750 -0.07709081  0.212090805 0.7240821
## IV-II   0.03625 -0.10834081  0.180840805 0.9733269
## V-II    0.01375 -0.13084081  0.158340805 0.9997181
## VI-II  -0.07250 -0.21709081  0.072090805 0.6611505
## IV-III -0.03125 -0.17584081  0.113340805 0.9861403
## V-III  -0.05375 -0.19834081  0.090840805 0.8705387
## VI-III -0.14000 -0.28459081  0.004590805 0.0624025
## V-IV   -0.02250 -0.16709081  0.122090805 0.9969635
## VI-IV  -0.10875 -0.25334081  0.035840805 0.2356038
## VI-V   -0.08625 -0.23084081  0.058340805 0.4817411

plot(eggs_posthoc)

Simplification d’une ANOVA à deux facteurs

# Simplification des données
eggs %>.%
  sgroup_by(., Technician) %>.%
  ssummarise(.,
    Fat_mean = fmean(Fat),
    Lab      = funique(Lab)) ->
  eggs_means
skimr::skim(eggs_means)

Data summary
Name	eggs_means
Number of rows	12
Number of columns	3
_______________________
Column type frequency:
factor	2
numeric	1
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
Technician	0	1	FALSE	12	I.o: 1, II.: 1, III: 1, IV.: 1
Lab	0	1	FALSE	6	I: 2, II: 2, III: 2, IV: 2

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Fat_mean	0	1	0.39	0.13	0.17	0.36	0.37	0.39	0.72	▁▇▂▁▁

# Résumé des données simplifiées
eggs_means %>.%
  sgroup_by(., Lab) %>.%
  ssummarise(.,
    moyenne  = fmean(Fat_mean),
    variance = fsd(Fat_mean),
    n        = fnobs(Fat_mean))

## # A data.trame: [6 × 4]
##   Lab   moyenne variance     n
##   <fct>   <dbl>    <dbl> <int>
## 1 I       0.58   0.202       2
## 2 II      0.34   0.0354      2
## 3 III     0.408  0.0530      2
## 4 IV      0.376  0.00177     2
## 5 V       0.354  0.00884     2
## 6 VI      0.267  0.131       2

# Représentation graphique adéquate de nos données résumées
chart(eggs_means, Fat_mean ~ Lab) +
  geom_point() +
  stat_summary(geom = "point", fun = "mean", color = "red", size = 3)

# Test de Bartlett pour l'homoscédasticité
bartlett.test(data = eggs_means, Fat_mean ~ Lab)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Fat_mean by Lab
## Bartlett's K-squared = 10.452, df = 5, p-value = 0.0634

# ANOVA à un facteur sur les données moyennées
eggs_means_anova <- lm(data = eggs_means, Fat_mean ~ Lab)
anova(eggs_means_anova) %>.%
  tabularise(.)

Analyse de la variance Réponse : Fat_mean
Terme	Ddl	Somme des carrés	Carrés moyens	Valeur de Fobs.	Valeur de p
Lab	5	0.1108	0.0222	2.15	0.19
Résidus	6	0.0619	0.0103
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05

# Graphique quantile-quantile des résidus de notre modèle simplifié
chart$qqplot(eggs_means_anova)

ANOVA avec effet aléatoire (split-plot)

# Lecture des données Penicillin depuis le package lme4
pen <- read("Penicillin", package = "lme4")
# Premières et dernières lignes du jeu de données
tabularise$headtail(pen, n = 16)

Diamètre de la zone d'inhibition [mm]	Boîte de Petri	Échantillon de pénicilline
27	a	A
23	a	B
26	a	C
23	a	D
23	a	E
21	a	F
27	b	A
23	b	B
...	...	...
24	w	E
19	w	F
24	x	A
21	x	B
24	x	C
22	x	D
21	x	E
18	x	F
Premières et dernières 8 lignes d'un total de 144

# Détermination du nombre de replicats
replications(data = pen, diameter ~ sample + plate)

## sample  plate 
##     24      6

# Visualisation des données
chart(data = pen, diameter ~ plate | sample) +
  geom_point()

# Tests de Batlett en fonction de sample et de plate
bartlett.test(data = pen, diameter ~ sample)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  diameter by sample
## Bartlett's K-squared = 0.56002, df = 5, p-value = 0.9898

bartlett.test(data = pen, diameter ~ plate)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  diameter by plate
## Bartlett's K-squared = 3.4151, df = 23, p-value = 1

# ANOVA classique à deux facteurs sans interactions
pen_anova <- lm(data = pen, diameter ~ sample + plate)
anova(pen_anova) %>.%
  tabularise(.)

Analyse de la variance Réponse : diameter
Terme	Ddl	Somme des carrés	Carrés moyens	Valeur de Fobs.	Valeur de p
sample	5	449.2	89.844	297.1	< 2·10-16	***
plate	23	105.9	4.604	15.2	< 2·10-16	***
Résidus	115	34.8	0.302
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05

# Graphique quantile-quantile des résidus du modèle
chart$qqplot(pen_anova)

# Graphique des résidus en fonction des valeurs prédites par le modèle
chart$resfitted(pen_anova)

# Ajustement du modèle en parcelles divisées avec lmeTest::lmer()
pen_split_plot <- lmerTest::lmer(data = pen, diameter ~ sample + (1 | plate))
pen_split_plot

## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: diameter ~ sample + (1 | plate)
##    Data: pen
## REML criterion at convergence: 308.2793
## Random effects:
##  Groups   Name        Std.Dev.
##  plate    (Intercept) 0.8467  
##  Residual             0.5499  
## Number of obs: 144, groups:  plate, 24
## Fixed Effects:
## (Intercept)      sampleB      sampleC      sampleD      sampleE      sampleF  
##      25.167       -3.208       -0.250       -2.292       -2.208       -5.208

# Tableau ANOVA de notre modèle split-plot
anova(pen_split_plot) %>.%
  tabularise(.)


Terme	Somme des carrés	Carrés moyens	Ddl. Num.	Ddl dénom.	Valeur de Fobs.	Valeur de p
sample	449	89.8	5	115	297	< 2·10-16	***
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05

# Intervalles de confiance à 95% sur les paramètres du modèle split-plot
confint(pen_split_plot)

## Computing profile confidence intervals ...

##                  2.5 %      97.5 %
## .sig01       0.6243175  1.15107599
## .sigma       0.4768022  0.61436044
## (Intercept) 24.7615793 25.57175404
## sampleB     -3.5153786 -2.90128804
## sampleC     -0.5570453  0.05704529
## sampleD     -2.5987120 -1.98462138
## sampleE     -2.5153786 -1.90128804
## sampleF     -5.5153786 -4.90128804

# Comparaisons multiples HSD de Tukey sur sample de notre modèle split-plot
summary(pen_posthoc <- confint(multcomp::glht(pen_split_plot,
  linfct = multcomp::mcp(sample = "Tukey"))))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmerTest::lmer(formula = diameter ~ sample + (1 | plate), data = pen)
## 
## Linear Hypotheses:
##            Estimate Std. Error z value Pr(>|z|)    
## B - A == 0 -3.20833    0.15875 -20.210   <1e-04 ***
## C - A == 0 -0.25000    0.15875  -1.575    0.615    
## D - A == 0 -2.29167    0.15875 -14.436   <1e-04 ***
## E - A == 0 -2.20833    0.15875 -13.911   <1e-04 ***
## F - A == 0 -5.20833    0.15875 -32.809   <1e-04 ***
## C - B == 0  2.95833    0.15875  18.635   <1e-04 ***
## D - B == 0  0.91667    0.15875   5.774   <1e-04 ***
## E - B == 0  1.00000    0.15875   6.299   <1e-04 ***
## F - B == 0 -2.00000    0.15875 -12.598   <1e-04 ***
## D - C == 0 -2.04167    0.15875 -12.861   <1e-04 ***
## E - C == 0 -1.95833    0.15875 -12.336   <1e-04 ***
## F - C == 0 -4.95833    0.15875 -31.234   <1e-04 ***
## E - D == 0  0.08333    0.15875   0.525    0.995    
## F - D == 0 -2.91667    0.15875 -18.373   <1e-04 ***
## F - E == 0 -3.00000    0.15875 -18.898   <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

oma <- par(oma = c(0, 5.1, 0, 0))
plot(pen_posthoc)

par(oma)
rm(oma)

# Graphique quantile-quantile du modèle split-plot
library(broom.mixed) # Required for mixed models
#chart$qqplot(pen_split_plot) # Not implemented yet
pen_split_plot %>.%
  augment(.) %>.%
  car::qqPlot(.$.resid, distribution = "norm",
    envelope = 0.95, col = "Black", xlab = "Quantiles théoriques",
    ylab = "Résidus standardisés")

## [1] 137  14

# Graphique des résidus en fonction des valeurs prédites par le modèle split-plot
#chart$resfitted(pen_split_plot) # Not implemented yet
pen_split_plot %>.%
  augment(.) %>.%
  chart(., .resid ~ .fitted) +
  geom_point() +
  geom_hline(yintercept = 0) +
  geom_smooth(se = FALSE, method = "loess", formula = y ~ x) +
  labs(x = "Valeurs prédites", y = "Résidus") +
  ggtitle("Distribution des résidus - pen_split_plot")

# Modèle avec deux facteurs aléatoires
pen_split_plot2 <- lmerTest::lmer(data = pen, diameter ~ (1 | sample) + (1 | plate))
pen_split_plot2

## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: diameter ~ (1 | sample) + (1 | plate)
##    Data: pen
## REML criterion at convergence: 330.8606
## Random effects:
##  Groups   Name        Std.Dev.
##  plate    (Intercept) 0.8467  
##  sample   (Intercept) 1.9316  
##  Residual             0.5499  
## Number of obs: 144, groups:  plate, 24; sample, 6
## Fixed Effects:
## (Intercept)  
##       22.97

# Intervalles de confiance pour notre modèle avec deux fateurs aléatoires
confint(pen_split_plot2)

## Computing profile confidence intervals ...

##                  2.5 %     97.5 %
## .sig01       0.6335660  1.1821040
## .sig02       1.0957894  3.5562912
## .sigma       0.4858454  0.6294535
## (Intercept) 21.2666274 24.6778176

ANOVA avec effet aléatoire (mesures répétées)

# Lecture des données sleep du package lme4
sleep <- read("sleepstudy", package = "lme4")
sleep

## # A data.trame: [180 × 3]
##    reaction  days subject
##       <dbl> <dbl> <fct>  
##  1     250.     0 308    
##  2     259.     1 308    
##  3     251.     2 308    
##  4     321.     3 308    
##  5     357.     4 308    
##  6     415.     5 308    
##  7     382.     6 308    
##  8     290.     7 308    
##  9     431.     8 308    
## 10     466.     9 308    
## # ℹ 170 more rows

# Exploration des données
skimr::skim(sleep)

Data summary
Name	sleep
Number of rows	180
Number of columns	3
_______________________
Column type frequency:
factor	1
numeric	2
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
subject	0	1	FALSE	18	308: 10, 309: 10, 310: 10, 330: 10

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
reaction	0	1	298.51	56.33	194.33	255.38	288.65	336.75	466.35	▃▇▆▂▁
days	0	1	4.50	2.88	0.00	2.00	4.50	7.00	9.00	▇▇▇▇▇

# Graphique adéquat pour visualiser des données pour modèles à mesures répétées
chart(data = sleep, reaction ~ days %col=% subject) +
  geom_line()

# Autre graphique adéquat avec des facettes
chart(data = sleep, reaction ~ days | subject) +
  geom_line()

# Graphique à facette avec représentation d'une tendance linéaire
chart(data = sleep, reaction ~ days | subject) +
  geom_point() +
  stat_smooth(method = "lm") # Ajuste une droite sur les données

## `geom_smooth()` using formula = 'y ~ x'

# Modèle à mesurées répétées avec lmerTest::lmer()
sleep_rep <- lmerTest::lmer(data = sleep, reaction ~ days + (days | subject))
sleep_rep

## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: reaction ~ days + (days | subject)
##    Data: sleep
## REML criterion at convergence: 1743.628
## Random effects:
##  Groups   Name        Std.Dev. Corr
##  subject  (Intercept) 24.741       
##           days         5.922   0.07
##  Residual             25.592       
## Number of obs: 180, groups:  subject, 18
## Fixed Effects:
## (Intercept)         days  
##      251.41        10.47

# Tableau de l'ANOVA de notre modèle à mesures répétées
anova(sleep_rep) %>.%
  tabularise(.)

Analyse de la variance de type III avec méthode Sattertwaite
Terme	Somme des carrés	Carrés moyens	Ddl. Num.	Ddl dénom.	Valeur de Fobs.	Valeur de p
days	30031	30031	1	17	45.9	3.26·10-06	***
0 <= '*' < 0.001 < '' < 0.01 < '*' < 0.05

# Résumé du modèle
summary(sleep_rep)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: reaction ~ days + (days | subject)
##    Data: sleep
## 
## REML criterion at convergence: 1743.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9536 -0.4634  0.0231  0.4634  5.1793 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  subject  (Intercept) 612.10   24.741       
##           days         35.07    5.922   0.07
##  Residual             654.94   25.592       
## Number of obs: 180, groups:  subject, 18
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  251.405      6.825  17.000  36.838  < 2e-16 ***
## days          10.467      1.546  17.000   6.771 3.26e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## days -0.138

# Intervalles de confiance à 95% sur les paramètres du modèle à mesures répétées
confint(sleep_rep)

## Computing profile confidence intervals ...

##                   2.5 %      97.5 %
## .sig01       14.3814778  37.7139771
## .sig02       -0.4815007   0.6849858
## .sig03        3.8011641   8.7533828
## .sigma       22.8982669  28.8579965
## (Intercept) 237.6806955 265.1295146
## days          7.3586533  13.5759187

# Graphique quantile-quantile des résidus du modèle à mesures répétées
library(broom.mixed) # Required for mixed models
#chart$qqplot(sleep_rep) # Not implemented yet
sleep_rep %>.%
  augment(.) %>.%
  car::qqPlot(.$.resid, distribution = "norm",
    envelope = 0.95, col = "Black", xlab = "Quantiles théoriques",
    ylab = "Résidus standardisés")

## [1] 57  8

# Distribution des résidus du modèle à mesures répétées par les valeurs prédites
#chart$resfitted(sleep_rep) # Not implemented yet
sleep_rep %>.%
  augment(.) %>.%
  chart(., .resid ~ .fitted) +
  geom_point() +
  geom_hline(yintercept = 0) +
  geom_smooth(se = FALSE, method = "loess", formula = y ~ x) +
  labs(x = "Valeurs prédites", y = "Résidus") +
  ggtitle("Distribution des résidus - pen_split_plot")

Syntaxe de R

# Un vecteur contenant 4 valeurs numériques nommées a, b, c et d
v <- c(a = 2.6, b = 7.1, c = 4.9, d = 5.0)
# Le second élément du vecteur
v[2]

##   b 
## 7.1

# Le premier et le troisième élément
v[c(1, 3)]

##   a   c 
## 2.6 4.9

# Les 3 premiers éléments avec la séquence 1, 2, 3 issue de 1:3
v[1:3]

##   a   b   c 
## 2.6 7.1 4.9

# Tout le vecteur v, sauf le 2ème élément
v[-2]

##   a   c   d 
## 2.6 4.9 5.0

# Élimination du 2ème et du 3ème élément
v[-(2:3)]

##   a   d 
## 2.6 5.0

# Élimination du dernier élément
v[-length(v)]

##   a   b   c 
## 2.6 7.1 4.9

# Élément de v s'appelant 'a'
v['a']

##   a 
## 2.6

# Les éléments 'b' et 'd'
v[c('b', 'd')]

##   b   d 
## 7.1 5.0

# Garder le premier et le quatrième élément
v[c(TRUE, FALSE, FALSE, TRUE)]

##   a   d 
## 2.6 5.0

# Garder le premier et le troisième (recyclage des indices une seconde fois)
v[c(TRUE, FALSE)]

##   a   c 
## 2.6 4.9

# Déterminer quel élément est plus grand que 3 dans v
v > 3

##     a     b     c     d 
## FALSE  TRUE  TRUE  TRUE

# Utilisation de cette instruction comme indiçage pour filtrer les éléments de v
v[v > 3]

##   b   c   d 
## 7.1 4.9 5.0

# Filtrage d'un vecteur en R de base
v[v > 3 & v <= 5]

##   c   d 
## 4.9 5.0

# Création d'un data frame
df <- dtx_rows(
  ~x, ~y, ~z,
   1,  2,  3,
   4,  5,  6
)
df

## # A data.trame: [2 × 3]
##       x     y     z
##   <dbl> <dbl> <dbl>
## 1     1     2     3
## 2     4     5     6

# Élément à la première ligne, colonnes 2 et 3
df[1, 2:3]

## # A data.trame: [1 × 2]
##       y     z
##   <dbl> <dbl>
## 1     2     3

# Toute la seconde ligne
df[2, ]

## # A data.trame: [1 × 3]
##       x     y     z
##   <dbl> <dbl> <dbl>
## 1     4     5     6

# Toute la seconde colonne
df[ , 2]

## # A data.trame: [2 × 1]
##       y
##   <dbl>
## 1     2
## 2     5

# Tout le tableau (pas très utile !)
df[ , ]

## # A data.trame: [2 × 3]
##       x     y     z
##   <dbl> <dbl> <dbl>
## 1     1     2     3
## 2     4     5     6

# Lignes pour lesquelles x est plus grand que 3 et colonnes nommée 'y' et 'z'
df[df$x > 3, c('y', 'z')]

## # A data.trame: [1 × 2]
##       y     z
##   <dbl> <dbl>
## 1     5     6

# Récupération d'une colonne d'un data frame
df[[2]]

## [1] 2 5

df[['y']]

## [1] 2 5

df$y

## [1] 2 5

# Remplacer la troisième colonne par des nouvelles valeurs
df[ , 3] <- c(-10, -15)
df

## # A data.trame: [2 × 3]
##       x     y     z
##   <dbl> <dbl> <dbl>
## 1     1     2   -10
## 2     4     5   -15

# Ceci donne le même résultat
df$z <- c(-10, -15)
df

## # A data.trame: [2 × 3]
##       x     y     z
##   <dbl> <dbl> <dbl>
## 1     1     2   -10
## 2     4     5   -15

# Dialecte SciViews::R pour lire un jeu de données cas par variables
SciViews::R
# Lecture des données zooplankton
zoo <- read("zooplankton", package = "data.io", lang = "FR")
# Exploration des données
skimr::skim(zoo)

Data summary
Name	zoo
Number of rows	1262
Number of columns	20
_______________________
Column type frequency:
factor	1
numeric	19
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
class	0	1	FALSE	17	Cal: 288, Poe: 158, Déc: 126, Mal: 121

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
ecd	1	0.81	0.51	0.28	0.54	0.67	0.88	7.25	▇▁▁▁▁
area	1	0.72	1.74	0.06	0.23	0.35	0.61	41.27	▇▁▁▁▁
perimeter	1	6.45	6.31	0.95	2.82	4.33	7.81	89.26	▇▁▁▁▁
feret	1	1.81	1.55	0.31	0.93	1.34	2.05	10.59	▇▁▁▁▁
major	1	1.33	1.28	0.29	0.74	0.91	1.38	10.01	▇▁▁▁▁
minor	1	0.56	0.37	0.18	0.36	0.48	0.63	5.91	▇▁▁▁▁
mean	1	0.21	0.12	0.02	0.13	0.20	0.27	0.78	▆▇▂▁▁
mode	1	0.07	0.20	0.00	0.01	0.02	0.03	1.02	▇▁▁▁▁
min	1	0.01	0.00	0.00	0.00	0.00	0.01	0.03	▇▁▁▁▁
max	1	0.70	0.22	0.07	0.54	0.68	0.89	1.02	▁▃▇▇▇
std_dev	1	0.17	0.07	0.01	0.12	0.17	0.22	0.42	▃▇▆▂▁
range	1	0.69	0.22	0.06	0.53	0.68	0.88	1.02	▁▃▇▇▇
size	1	0.94	0.73	0.28	0.59	0.71	1.02	7.40	▇▁▁▁▁
aspect	1	0.54	0.24	0.06	0.35	0.52	0.74	1.00	▃▇▇▆▅
elongation	1	17.44	17.17	1.00	5.19	11.93	24.25	134.36	▇▂▁▁▁
compactness	1	6.24	5.43	1.10	2.35	4.46	8.37	43.41	▇▂▁▁▁
transparency	1	0.09	0.11	0.00	0.01	0.05	0.12	0.54	▇▂▁▁▁
circularity	1	0.31	0.24	0.02	0.12	0.22	0.43	0.91	▇▅▂▁▂
density	1	0.13	0.25	0.00	0.04	0.07	0.15	5.74	▇▁▁▁▁

# Tableau de contingence en R de base
table(zoo$class)

## 
##          Annélide    Appendiculaire         Calanoïde      Chaetognathe         Cirripède         Cladocère          Cnidaire        Cyclopoïde 
##                50                36               288                51                22                50                22                50 
##          Décapode      Oeuf_allongé         Oeuf_rond           Poisson       Gastéropode     Harpacticoïde      Malacostracé Poecilostomatoïde 
##               126                50                49                50                50                39               121               158 
##          Protiste 
##                50

# Tbleau de contingence avec interface formule
mosaic::tally(data = zoo, ~ class)

## class
##          Annélide    Appendiculaire         Calanoïde      Chaetognathe         Cirripède         Cladocère          Cnidaire        Cyclopoïde 
##                50                36               288                51                22                50                22                50 
##          Décapode      Oeuf_allongé         Oeuf_rond           Poisson       Gastéropode     Harpacticoïde      Malacostracé Poecilostomatoïde 
##               126                50                49                50                50                39               121               158 
##          Protiste 
##                50

# Tableau de contingence en syntaxe Tidyverse
zoo %>%
  group_by(class) %>%
  summarise(n())

## # A tibble: 17 × 2
##    class             `n()`
##    <fct>             <int>
##  1 Annélide             50
##  2 Appendiculaire       36
##  3 Calanoïde           288
##  4 Chaetognathe         51
##  5 Cirripède            22
##  6 Cladocère            50
##  7 Cnidaire             22
##  8 Cyclopoïde           50
##  9 Décapode            126
## 10 Oeuf_allongé         50
## 11 Oeuf_rond            49
## 12 Poisson              50
## 13 Gastéropode          50
## 14 Harpacticoïde        39
## 15 Malacostracé        121
## 16 Poecilostomatoïde   158
## 17 Protiste             50

# Fonction Tidyverse spécialisée pour le contingentement
count(zoo, class)

## # A data.trame: [17 × 2]
##    class                 n
##    <fct>             <int>
##  1 Annélide             50
##  2 Appendiculaire       36
##  3 Calanoïde           288
##  4 Chaetognathe         51
##  5 Cirripède            22
##  6 Cladocère            50
##  7 Cnidaire             22
##  8 Cyclopoïde           50
##  9 Décapode            126
## 10 Oeuf_allongé         50
## 11 Oeuf_rond            49
## 12 Poisson              50
## 13 Gastéropode          50
## 14 Harpacticoïde        39
## 15 Malacostracé        121
## 16 Poecilostomatoïde   158
## 17 Protiste             50

# Filtrage de lignes et sélection de colonnes avec [,] en R de base
zoo2 <- zoo[zoo$class == "Oeuf_allongé" | zoo$class == "Oeuf_rond",
  c("aspect", "area", "class")]
zoo2

## # A data.trame: [99 × 3]
##    aspect  area class       
##     <dbl> <dbl> <fct>       
##  1  0.965 0.385 Oeuf_rond   
##  2  0.976 6.04  Oeuf_rond   
##  3  0.964 0.271 Oeuf_rond   
##  4  0.911 0.259 Oeuf_rond   
##  5  0.950 0.322 Oeuf_rond   
##  6  0.959 0.335 Oeuf_rond   
##  7  0.976 0.343 Oeuf_rond   
##  8  0.975 1.11  Oeuf_rond   
##  9  0.449 0.495 Oeuf_allongé
## 10  0.987 2.02  Oeuf_rond   
## # ℹ 89 more rows

# Filtrage de lignes et sélection de colonnes avec le Tidyverse
zoo %>.%
  dplyr::filter(., class == "Oeuf_allongé" | class == "Oeuf_rond") %>.%
  dplyr::select(., aspect, area, class) ->
  zoo2
zoo2

## # A data.trame: [99 × 3]
##    aspect  area class       
##     <dbl> <dbl> <fct>       
##  1  0.965 0.385 Oeuf_rond   
##  2  0.976 6.04  Oeuf_rond   
##  3  0.964 0.271 Oeuf_rond   
##  4  0.911 0.259 Oeuf_rond   
##  5  0.950 0.322 Oeuf_rond   
##  6  0.959 0.335 Oeuf_rond   
##  7  0.976 0.343 Oeuf_rond   
##  8  0.975 1.11  Oeuf_rond   
##  9  0.449 0.495 Oeuf_allongé
## 10  0.987 2.02  Oeuf_rond   
## # ℹ 89 more rows

# Calcul d'une nouvelle variable en R de base
zoo2$log_area <- log10(zoo2$area)
head(zoo2)

## # A data.trame: [6 × 4]
##   aspect  area class     log_area
##    <dbl> <dbl> <fct>        <dbl>
## 1  0.965 0.385 Oeuf_rond   -0.414
## 2  0.976 6.04  Oeuf_rond    0.781
## 3  0.964 0.271 Oeuf_rond   -0.568
## 4  0.911 0.259 Oeuf_rond   -0.587
## 5  0.950 0.322 Oeuf_rond   -0.492
## 6  0.959 0.335 Oeuf_rond   -0.475

# Calcul d'une nouvelle variable avec mutate() du Tidyverse
zoo2 <- mutate(zoo2, log_area = log10(area))
head(zoo2)

## # A data.trame: [6 × 4]
##   aspect  area class     log_area
##    <dbl> <dbl> <fct>        <dbl>
## 1  0.965 0.385 Oeuf_rond   -0.414
## 2  0.976 6.04  Oeuf_rond    0.781
## 3  0.964 0.271 Oeuf_rond   -0.568
## 4  0.911 0.259 Oeuf_rond   -0.587
## 5  0.950 0.322 Oeuf_rond   -0.492
## 6  0.959 0.335 Oeuf_rond   -0.475

zoo2 <- mutate_(zoo2, log_area = ~log10(area))
head(zoo2)

## # A data.trame: [6 × 4]
##   aspect  area class     log_area
##    <dbl> <dbl> <fct>        <dbl>
## 1  0.965 0.385 Oeuf_rond   -0.414
## 2  0.976 6.04  Oeuf_rond    0.781
## 3  0.964 0.271 Oeuf_rond   -0.568
## 4  0.911 0.259 Oeuf_rond   -0.587
## 5  0.950 0.322 Oeuf_rond   -0.492
## 6  0.959 0.335 Oeuf_rond   -0.475

# Tous les niveaux sont toujours là
levels(zoo2$class)

##  [1] "Annélide"          "Appendiculaire"    "Calanoïde"         "Chaetognathe"      "Cirripède"         "Cladocère"         "Cnidaire"         
##  [8] "Cyclopoïde"        "Décapode"          "Oeuf_allongé"      "Oeuf_rond"         "Poisson"           "Gastéropode"       "Harpacticoïde"    
## [15] "Malacostracé"      "Poecilostomatoïde" "Protiste"

# Ne retenir que les niveaux relatifs aux œufs
zoo2$class <- droplevels(zoo2$class)
# C'est mieux !
levels(zoo2$class)

## [1] "Oeuf_allongé" "Oeuf_rond"

# Graphique en R de base
plot(zoo2$log_area, zoo2$aspect, col = zoo2$class)
legend("bottomright", legend = c("Oeuf allongé", "Oeuf rond"), col = 1:2, pch = 1)

# Idem, mais avec interface formule
plot(data = zoo2, aspect ~ log_area, col = class)
legend("bottomright", legend = c("Oeuf allongé", "Oeuf rond"), col = 1:2, pch = 1)

# Graphique lattice (xyplot) à partir de chart()
chart$xyplot(data = zoo2, aspect ~ log_area, groups = zoo2$class, auto.key = TRUE)

# Graphique ggplot2 (Tidyverse)
ggplot(data = zoo2, aes(x = log_area, y = aspect, col = class)) +
  geom_point()

# chart() avec aes() 
chart(data = zoo2, aes(x = log_area, y = aspect, col = class)) +
  geom_point()

# chart() avec une formule élargie
chart(data = zoo2, aspect ~ log_area %col=% class) +
  geom_point()

SDD I module 10 : Variance II

Ph. Grosjean et G. Engels

2025-2026

ANOVA à deux facteurs sans interactions

ANOVA à deux facteurs croisés complet

ANOVA à deux facteurs hiérarchisés

Simplification d’une ANOVA à deux facteurs

ANOVA avec effet aléatoire (split-plot)

ANOVA avec effet aléatoire (mesures répétées)

Syntaxe de R