diapos_seance3.Rmd

---
title: "ANOVA à 1 facteur"
subtitle: "Séance 3 de *modèles linéaires*"
author: "Florent Chuffart & Magali Richard"
date: "`r Sys.Date()`"
output: 
  slidy_presentation:

# output: html_document
---

```{r setup, include=FALSE}
options(htmltools.dir.version = FALSE, width = 75)
knitr::opts_chunk$set(echo = TRUE, fig.width=15, fig.heigh=5, fig.align='center', dev='png', dpi = 95,
echo=FALSE, results="hide")

```

---

# Organisation

## Pré-requis
 
 - R https://cran.r-project.org 
 - RStudio https://www.rstudio.com

## Evaluation

 - individuelle sur un jeu de données aprés 4 séances de cours 
 - data challenge à la fin des séances

## Supports

 - https://github.com/magrichard/cours_modlin_M1


## Découpage prévisionnelle

Séance 1 - 3h (10/01) Régression linéaire simple

Séance 2 - 3h (24/01) Régression linéaire multiple

Séance 3 - 3h (31/01) Anova 1 facteur / partie A

Séance 4 - 3h (14/02) Anova 1 facteur / partie B

Séance 5 - 3h (21/02) Anova 2 facteurs

Séance 6 - 3h (07/03) Anova multiple

Séance 7 - 3h (11/03) Anova facteurs emboités

Séance 8 - 3h (14/03) Ancova


---

# Plan

### ANOVA à 1 facteur

1. Rappels regression linéaire
2. Présentation de l’ANOVA à un facteur
3. Décomposition de la variance
4. Test de l’ANOVA
5. Comparaisons multiples

---

# Rappels regression linéaire

### Regression linéaire simple

**Objectif de la regression linéaire simple** Considérons $x$ et $Y$ deux variables quantitatives. Y est expliquée (modélisée) par  la variable explicative $x$. 

Modèle : $$\widehat{Y} = E(Y)  = f(x) = \beta_0 + \beta_1 x$$

$$Y_i = \beta_0 + \beta_1 x + e_i$$

```{r}
d = read.table("data/data_nutri.csv", header=TRUE, sep=",", row.names = 1)
layout(matrix(1:2, 1, byrow=TRUE), respect=TRUE)
plot(d$taille, d$poids, main="poids~taille")
## Model
# Y~X
# E(Y) = b.X
# E(Y) = b_0 + b_1.X
# Y_i = b_0 + b_1.X_i + e_i
m = lm(d$poids~d$taille)
# m$coefficients
# residuals
# m$residuals
t = jitter(d$taille)
suppressWarnings(arrows(t, d$poids, t, d$poids-m$residuals, col=adjustcolor(4, alpha.f=0.2), length=0.1))
suppressWarnings(arrows(t, d$poids-m$residuals, t, mean(d$poids), col=adjustcolor(2, alpha.f=0.2), length=0.1))
abline(a=m$coefficients[[1]], b=m$coefficients[[2]]) # /!\ y = b.x + a
abline(h=mean(d$poids), lty=2)
legend("topleft",c("regression line", "observations mean", "residual variance", "explained variance"), col=c(1,1,4,2), lty=c(1,2,1,1), cex=.8)
```

la méthode des *moindres carrés ordinaires (MCO)* minimise l’*erreur quadratique moyenne* (EQM, la somme du carré des résidus).

$$EQM(\beta_0, \beta_1) = \sum^{n}_{i=1} (Y_i - f(x_i))^2$$

La variance totale se décompose en la somme de la variance expliquée (par la regression) et la variance résiduelle.

$$ \sum(Y_i - \overline{Y}_n)^2 = \sum(\widehat{Y}_i - \overline{Y}_n)^2 + \sum(Y_i - \widehat{Y}_i)^2 $$

avec $Y_i$ la *i*-éme valeur observée pour $i=1,...,n$, $\widehat{Y}_i = f(x_i)$ la valeur estimée par le modèle pour $i=1,...,n$ et $\overline{Y}_n$  la moyenne empirique des $n$ observations.

$$ SC_{tot} = SC_{reg} + SC_{res}$$ 

Avec $SC_{tot}$ la somme des carrés totale,  $SC_{reg}$ la somme des carrés due à la régression et  $SC_{res}$ la somme des carrés des résidus.

Le coefficient de determination $R^2$ mesure du pourcentage de la variance expliquée par le modèle.

$$ R^2 = \frac{SC_{reg}}{SC_{tot}}$$


```{r echo=TRUE, results="verbatim"}
sc_tot = sum((d$poids - mean(d$poids))^2)
sc_reg = sum((d$poids-m$residuals - mean(d$poids))^2)
sc_res = sum(m$residuals^2)

sc_tot
sc_reg
sc_res

sc_reg / sc_tot

m = lm(d$poids~d$taille)
summary(m)
```

### Regression linéaire multiple

**Objectif de la regression linéaire multiple** Considérons $X_1, ..., X_p$ et $Y$ plusieurs variables quantitatives. Y est expliquée (modélisée) par  les variable explicative $X_1, ..., X_p$. 

Modèle : 

$$\widehat{Y}  = \beta_0 + \beta_1X_1 + ... + \beta_pX_p + \epsilon  $$


Utiliser les **moindres carrés** revient à minimiser la quantité suivante: 

$$min_{\beta_0, ...,\beta_p}\sum^n_{i=1}\Big(y_i -(\beta_0+ \beta_1X_1 + ... + \beta_pX_p)\Big)^2$$

---
 
# Présentation de l'ANOVA à un facteur

$$Y \sim X$$

*X* - La régression linéaire se caractérise par des variables explicatives **quantitatives** (Ex : taille, poids). 
L’analyse de la variance (ANOVA) se caractérise par des variables explicatives **qualitatives** (Ex : sexe, couleur des yeux).

*Y* - Dans les deux cas la variable expliquée est quantitatives.

L’**ANOVA** à un facteur est un à la fois un **modèle** statistique fondé sur la décomposition de la variance et **test** statistique permettant de comparer les moyennes de plusieurs variables aléatoires indépendantes, gaussiennes et de même variance.

Ex : le poids moyen de différent groupe d’individus.

L’analyse de la variance est l’une des procédures les plus utilisées dans les applications de la statistique ainsi que dans les méthodes d’analyse de données.






### Exemple *InsectSprays*

Cette base de données comprend le comptage du nombre d’abeilles présentent dans parcelle agricole (*count*, *numeric*) en fonction du type de spray utilisé sur cette parcelles (*spray*, *factor*).

```{r results="verbatim", echo=TRUE}
data("InsectSprays")
d = InsectSprays
head(d)
tail(d)
```

$$Y \sim X$$

On cherche à expliquer le nombre d’abeilles (Y) en fonction du type de spray (X).

On s’intérroge sur l’effet du type de spray utilisé (X) sur le nombre d’abeilles observé (Y).


```{r, echo=TRUE}
layout(matrix(1:2, 1, byrow=TRUE), respect=TRUE)
p = plot(d$spray, d$count, main="count~spray", xlab="spray", ylab="count", border="grey")
points(jitter(as.numeric(d$spray)), d$count)
```


### L'ANOVA à un facteur, un modèle statistique

$$Y_{ij} = \mu + \beta_i + \epsilon_{ij}$$


```{r}
y_i = sapply(levels(d$spray), function(s) {
 mean(d[d$spray==s,]$count)
})
layout(matrix(1:2, 1, byrow=TRUE), respect=TRUE)
sp = jitter(as.numeric(d$spray))
plot(sp, d$count, main="count~spray", las=2, type="p", xaxt="n", xlab="spray", ylab="count")
suppressWarnings(arrows(sp, d$count, sp, rep(y_i,each=12), col=adjustcolor(4, alpha.f=0.5), length=0.05))
points(y_i, pch=5, col=2)
suppressWarnings(arrows(1:6, mean(y_i), 1:6, y_i, col=adjustcolor(2, alpha.f=0.9), length=0.05))
abline(h=mean(y_i), lty=2)

axis(1, 1:length(levels(d$spray)), labels=levels(d$spray))
legend("top", c("Y_ij", "beta_i", "residuals", "mu"), pch=c(1,5,5,NA), col=c(1,2,4,1), lty=c(0,0,0,2), border=0)
```

où : $$\sum_{i=1}^I \beta_i = 0$$ 

Ainsi : $$\mu_i = \mu + \beta_i, i=1,...,I$$






### L'ANOVA à un facteur, un test de comparaison des moyennes

Sous certaines hypothèses, l’ANOVA devient un test de comparaison des moyennes des facteurs.

Nous nous proposons de tester l’hypothèse nulle : $$H_0, \mu_1 =\mu_2 =···=\mu_I$$

contre l’hypothèse alternative : $$H_1, \text{les moyennes } \mu_i \text{ ne sont pas toutes égales.}$$

La méthode statistique qui permet d’effectuer ce test est appelée l’analyse de la variance à un facteur (*one way analysis of variance*).


*Intuition 1* - Sous l’hypothèse nulle, et sous les hypothèses du test, l’interval de confiance (au risque $\alpha$) de la moyenne de chaque groupe doit contenir la moyenne de toutes les observations. Si ce n’est pas le cas on rejette $H_0$ et on accepte $H_1$ au risque $\alpha$. Si c’est le cas on conserve $H_0$ au risque $\beta$ à calculer. 

*Intuition 2* - Si la variable observée pour chaque groupe est normalement distribuée, l’intervalle de confiance à $95\%$ de la moyenne de chaque groupe est très proche de la moyenne de chaque groupe +/- 2 fois l’écart type de chaque groupe. 

Ici, ce n’est  pas le cas. Après vérification des hypothèse du test on pourrait tenter de rejetter $H_0$ et accepter $H_1$ au risque $\alpha$ de $5\%$.


```{r results="verbatim",}
y_i = sapply(levels(d$spray), function(s) {
 mean(d[d$spray==s,]$count)
})
s2_i = sapply(levels(d$spray), function(s) {
 var(d[d$spray==s,]$count)
})
s_i = sqrt(s2_i)
layout(matrix(1:2, 1, byrow=TRUE), respect=TRUE)
sp = jitter(as.numeric(d$spray))
plot(sp, d$count, main="count~spray", las=2, type="p", xaxt="n", xlab="spray", ylab="count")
abline(h=mean(y_i), lty=2)
sp = 1:length(levels(d$spray))
suppressWarnings(arrows(sp, y_i, sp, y_i + 2 * sqrt(s2_i), col=4, length=0.05, angle=135))
suppressWarnings(arrows(sp, y_i, sp, y_i - 2 * sqrt(s2_i), col=4, length=0.05, angle=135))
points(y_i, pch=3, col=2)
axis(1, 1:length(levels(d$spray)), labels=levels(d$spray))
legend("top", c("observations", "means per group", "means +/- 2sd", "global mean"), pch=c(1,3,4,NA), col=c(1,2,4,1), lty=c(0,0,0,2))
```


### Hypothèses du test

Les résidus $\widehat{e_{ij}}$ sont associés, sans en être des réalisations, aux variables erreurs $\epsilon_{ij}$ qui sont inobservables et satisfont aux 3 conditions suivantes :

1. Elles sont indépendantes
3. Elles sont de loi gaussienne (normalité)
2. Elles ont même variance $\sigma^2$ inconnue (homogénéité des variances ou homoscédasticité)


```{r}
y_i = sapply(levels(d$spray), function(s) {
 mean(d[d$spray==s,]$count)
})
layout(matrix(1:2, 1, byrow=TRUE), respect=TRUE)
sp = jitter(as.numeric(d$spray))
plot(sp, d$count, main="count~spray", las=2, type="p", xaxt="n", xlab="spray", ylab="count")
suppressWarnings(arrows(sp, d$count, sp, rep(y_i,each=12), col=adjustcolor(4, alpha.f=0.5), length=0.05))
points(y_i, pch=3, col=2)
# suppressWarnings(arrows(1:6, mean(y_i), 1:6, y_i, col=adjustcolor(2, alpha.f=0.9), length=0.05))
# abline(h=mean(y_i), lty=2)

axis(1, 1:length(levels(d$spray)), labels=levels(d$spray))
legend("top", c("observations", "means per group", "residuals"), pch=c(1,3,5), col=c(1,2,4))
```








Ce test se fonde sur la décomposition de la variance.








---

# Décomposition de la variance

### Notations

Pour les observations, nous utilisons deux indices : 

- $i$ indifie le *i*-éme type de spray avec $i=1,...,I$.
- $j$ indifie la *j*-éme observation (réplicat) pour chaque type de spray avec $j=1,...,J(i)$.

Les observations sont notées $Y_{ij}$.

Lorsque les groupes échantillons sont de même taille, nous disons que l’expérience est équilibrée.

```{r results="verbatim", echo=TRUE}
table(d$spray)
```

Si l’expérience est équilibrée $J(i) = J , \forall i$.

Ainsi, les observations sont notées $Y_{ij}$

La moyenne observée du groupe $i$ avec $i = 1,...,I$ s’écrit : 
$$\overline{Y}_i = \frac{1}{J}\sum^J_{j=1}Y_{ij}$$

La variance observée du groupe $i$ avec $i = 1,...,I$ s’écrit : 
$$s_i^2(Y) = \frac{1}{J}\sum^J_{j=1}(Y_{ij} -\overline{Y}_i)^2$$


*Remarque* : Cette dernière formule exprime la variance non corrigée. Très souvent, dans les ouvrages ou les logiciels, c’est la variance corrigée qui est utilisée : au lieu d’être divisée par J, la somme est divisée par J − 1.

Dans notre exemple, la moyenne, la variance et l’écart type de chaque groupe sont : 

```{r results="verbatim"}
y_i = sapply(levels(d$spray), function(s) {
 mean(d[d$spray==s,]$count)
})
s2_i = sapply(levels(d$spray), function(s) {
 var(d[d$spray==s,]$count)
})
s_i = sqrt(s2_i)
print(y_i)
print(s2_i)
print(s_i)
```


```{r results="verbatim",}
y_i = sapply(levels(d$spray), function(s) {
 mean(d[d$spray==s,]$count)
})
s2_i = sapply(levels(d$spray), function(s) {
 var(d[d$spray==s,]$count)
})
s_i = sqrt(s2_i)
layout(matrix(1:2, 1, byrow=TRUE), respect=TRUE)
sp = jitter(as.numeric(d$spray))
plot(sp, d$count, main="count~spray", las=2, type="p", xaxt="n", xlab="spray", ylab="count")
points(y_i, pch=4, col=2)
sp = 1:length(levels(d$spray))
suppressWarnings(arrows(sp, y_i, sp, y_i + 2 * sqrt(s2_i), col=adjustcolor(4, alpha.f=0.5), length=0.05, angle=135))
suppressWarnings(arrows(sp, y_i, sp, y_i - 2 * sqrt(s2_i), col=adjustcolor(4, alpha.f=0.5), length=0.05, angle=135))
axis(1, 1:length(levels(d$spray)), labels=levels(d$spray))
legend("bottomright", c("Y_ij", "Y_i", "Y_i +/- 2 s_i"), pch=c(1,4,4), col=c(1,2,4))
```

### Propriétés fondamentales

Le test est fondé sur deux propriétés des moyennes et des variances.

(1) La moyenne de toutes les observations est la moyenne des moyennes de chaque échantillon. Ceci s’écrit :

$$\overline{y} 
 = \frac{1}{n} \sum_{j=1}^J \sum_{i=1}^I y_{ij} 
 = \frac{1}{n} \sum_{i=1}^I \sum_{j=1}^J y_{ij}
 = \frac{1}{I} \sum_{i=1}^I \overline{y_i} $$


(2) La variance de toutes les observations est la somme de la variance des moyennes et de la moyenne des variances. Ceci s’écrit :

$$s^2(y) =  
 = \frac{1}{n} \sum_{j=1}^J \sum_{i=1}^I (y_{ij} - \overline{y})^2 
 = \frac{1}{I} \sum_{i=1}^I (\overline{y_i} - \overline{y})^2 + \frac{1}{I}\sum_{i=1}^I s_i^2(y) $$

On conserve la proriété vu dans la regression linéaire
 
$$ SC_{tot} = SC_{F} + SC_{res}$$ 

Avec $SC_{tot}$ la somme des carrés totale,  $SC_{F}$ la somme des carrés des facteurs et  $SC_{res}$ la somme des carrés des résidus.


```{r}
y_i = sapply(levels(d$spray), function(s) {
 mean(d[d$spray==s,]$count)
})
layout(matrix(1:2, 1, byrow=TRUE), respect=TRUE)
sp = jitter(as.numeric(d$spray))
plot(sp, d$count, main="count~spray", las=2, type="p", xaxt="n", xlab="spray", ylab="count")
suppressWarnings(arrows(sp, d$count, sp, rep(y_i,each=12), col=adjustcolor(4, alpha.f=0.5), length=0.05))
points(y_i, pch=5, col=2)
suppressWarnings(arrows(1:6, mean(y_i), 1:6, y_i, col=adjustcolor(2, alpha.f=0.9), length=0.05))
abline(h=mean(y_i), lty=2)

axis(1, 1:length(levels(d$spray)), labels=levels(d$spray))
legend("bottomright", c("Y_ij", "Y_i", "residuals"), pch=c(1,5,5), col=c(1,2,4))
```

Le coefficient de determination $R^2$ mesure du pourcentage de la variance expliquée par le modèle.

$$ R^2 = \frac{SC_{F}}{SC_{tot}}$$

```{r echo=TRUE, results="verbatim"}
m = lm(d$count~d$spray)
summary(m)
```

### En pratique...

  - Réaliser plusieurs modélisations du jeu de données CO2.
  - Evaluer la puissance des modèles que vous proposez.

https://github.com/fchuffar/tp_data_co2

---

# Test de l'ANOVA

*Intuition* - Si l’hypothèse nulle $H_0$ est vraie alors la quantité $SC_F$ doit être petite par rapport à la quantité $SC_{res}$.
Par contre, si l’hypothèse alternative $H_1$ est vraie alors la quantité $SC_F$ doit être grande par rapport à la quantité $SC_{res}$.

```{r}
y_i = sapply(levels(d$spray), function(s) {
 mean(d[d$spray==s,]$count)
})
layout(matrix(1:2, 1, byrow=TRUE), respect=TRUE)
sp = jitter(as.numeric(d$spray))
plot(sp, d$count, main="count~spray", las=2, type="p", xaxt="n", xlab="spray", ylab="count")
suppressWarnings(arrows(sp, d$count, sp, rep(y_i,each=12), col=adjustcolor(4, alpha.f=0.5), length=0.05))
points(y_i, pch=5, col=2)
suppressWarnings(arrows(1:6, mean(y_i), 1:6, y_i, col=adjustcolor(2, alpha.f=0.9), length=0.05))
abline(h=mean(y_i), lty=2)

axis(1, 1:length(levels(d$spray)), labels=levels(d$spray))
legend("bottomright", c("Y_ij", "SC_F", "SC_res"), pch=c(1,5,5), col=c(1,2,4))
```

Pour comparer ces quantités, R. A. Fisher, après les avoir
"corrigées" par leurs degrés de liberté (ddl), a considéré leur rapport.

Nous appelons *carré moyen associé au facteur* le terme 

$$CM_F = \frac{SC_F}{I-1}$$

et *carré moyen résiduel* le terme :

$$CM_{res} = \frac{SC_{res}}{n-I}$$

Le carré moyen résiduel est un estimateur sans biais de la variance des erreurs  $\sigma^2$.
C’est pourquoi il est souvent également appelé variance résiduelle et presque systématiquement noté $S_{res}^2$ lorsqu’il sert à estimer la variance des erreurs.
Sa valeur observée sur l’échantillon est ainsi notée $cm_{res}$ ou $s_{res}^2$ .


Si les trois conditions sont satisfaites et si l’hypothèse nulle
$H_0$ est vraie alors

$$ F_{obs} = \frac{cm_F}{cm_{res}}$$

est une réalisation d’une variable aléatoire F qui suit une loi de Fisher à $I-1$ degrés de liberté au numérateur et $n-I$ degrés de liberté au dénominateur. Cette loi est notée  $\mathcal{F}_{I-1,n-I}$.

### Tableau de l’analyse de la variance 


```{r echo=TRUE, results="verbatim"}
m = lm(d$count~d$spray)
anova(m)
```

| Source de variation   | sc        | ddl  |  cm  | $F_{obs}$
| :-------------------- | :-------: | :--: | :--: | ---:
| Due au facteur        | $sc_{F}$  | $I-1$  | $cm_{F}=\frac{sc_{F}}{I-1}$ | $\frac{cm_{F}}{cm_{res}}$
| Résiduelle            | $sc_{res}$| $n-I$  | $cm_{res}=\frac{sc_{res}}{n-I}$ 
| Totale                | $sc_{tot}$| $n-1$ 










### Hypothèses du test

#### Indépendance


Il n’existe pas, dans un contexte général, de test statistique simple permettant d’étudier l’indépendance.
Ce sont les conditions de l’expérience qui nous permettront d’affirmer que nous sommes dans le cas de l’indépendance.


#### Normalité

Nous ne pouvons pas, en général, la tester pour chaque échantillon. En effet le nombre d’observations est souvent très limité pour chaque échantillon.
Nous allons donc la tester sur l’ensemble des données.

Remarquons que si les conditions sont satisfaites et si nous notons :
$$\epsilon_{ij} = Y_{ij} - \mu_i$$
 
alors c’est la même loi pour l’ensemble des unités.
$$ L(\epsilon_{ij}) = N(0 ; \sigma_2) $$

Les moyennes $\mu_i$ étant inconnues, nous les estimons par les
estimateurs de la moyenne : les $Y_i$.

Nous pouvons alors tester la normalité, avec le test de Shapiro-Wilk ou avec le test de Shapiro-Francia sur l’ensemble des résidus.



#### Homoscedasticité

Plusieurs tests permettent de tester l’égalité de plusieurs variances. Parmi ceux-ci, un test souvent utilisé est le test de Bartlett dont le protocole est le suivant.

L’hypothèse nulle : $$H_0, \sigma_1^2 = \sigma_2^2 = ...= \sigma_I^2$$

contre l’hypothèse alternative : $$H_1, \text{Les variances } \sigma_i^2 \text{ ne sont pas toutes égales.}$$


En pratique, nous pouvons l’appliquer lorsque les effectifs $n_i$ des $I$ échantillons sont tous au moins égaux à 3.

Ce test requiert la normalité des erreurs.

Note : nous avons également testé l’homoscedasticité avec *Goldfeld-Quandt* dans le cadre du premier TP et lors de la seconde séance.



### En pratique...

  - Mettre en place une analyse de variance pour étudier l’e􏰒ffet du régime sur les temps de coagulation.
  - Analyse graphique des donnée
  - Vérifier les hypothèses du test. 

```{r echo=TRUE, result="verbatim"}
# Loading data
d = faraway::coagulation
head(d)
table(d$diet)

# Graphicals
layout(matrix(1:2, 1, byrow=TRUE), respect=TRUE)
plot(d$diet, d$coag, main="coag~diet", xlab="diet", ylab="coag", border="grey")
points(jitter(as.numeric(d$diet)), d$coag)

# Test hypothesis
m = lm(coag~diet, d)
m2 = aov(coag~diet, d)

## normality
m$residuals
m2$residuals
m$residuals == m2$residuals
shapiro.test(m$residuals)
# H_0 normality
# H_1 not normality
# pval = 0.86 > 0.05, we keep H_0, residuals are normal

## homoscedasticity
bartlett.test(coag~diet, d)
# H_0 homoscedasticity
# H_1 not homoscedasticity
# pval = 0.6441 > 0.05, we keep H_0, variance of groups are equal

anova(m)

## ANOVA
# H_0 means of each group are equals
# H_1 means of each group  are not equals
# alpha = 5%
# pval = 4.6e-05 < 5% we reject H_0 and accept H_1 at threshold of 5%, means are not equals.


sc_F = 228
sc_res = 112

mc_F = sc_F/3
mc_F
mc_res = sc_res/20
mc_res

F_obs = mc_F/mc_res
F_obs

1 - pf(13.571, 3, 20)



```

---

# Comparaisons multiples


Lorsque que nous rejettons $H_0$, nous pouvons chercher à analyser les différences entre les groupes. Nous procédons alors à des tests qui vont répondre à la question suivante : D’où vient la différence ? Quelles moyennes sont différentes ?

Ces tests qui vont répondre à cette question sont les tests de comparaisons multiples, des adaptations du test de Student.

La plus connues est la correction de *Bonferroni*. Cela consiste à diviser le seuil par le nombre de comparaisons. Bonferroni a montré que cette procédure garantit un taux d’erreur global plus faible que le seuil initial.  

Dans notre exemple, nous avons décidé que les moyennes théoriques sont différentes dans leur ensemble, mais nous aurions très bien pu trouver le contraire.
Comme nous avons décidé que les moyennes théoriques sont différentes dans leur ensemble que le facteur étudié est à effets fixes et qu’il a plus de trois modalités, nous pourrions essayer de déterminer là où résident les différences avec un des tests de comparaisons multiples.

```{r echo=TRUE, results="verbatim"}
summary(m)
```

---

# Réferences

Ce cours s’inspire des références suivantes :

- Frédéric Bertrand & Myriam Maumy-Bertrand
- Franck Picard
- Irenne Ganaz

http://irma.math.unistra.fr/~fbertran/enseignement/Master1_2016/Master1_Cours1.pdf