Analyse de la variance à 2 facteurs

# Analyse de la variance à 2 facteurs
### Marie-Pierre Etienne
### <a href="mailto:marie-pierre.etienne@agrocampus-ouest.fr" class="email">marie-pierre.etienne@agrocampus-ouest.fr</a>
### 2021/11/19 (updated: 2023-09-08)

---

# Introduction

## Etude de la fréquence cardiaque au repos

On a mesuré la fréquence cardiaque de 12 femmes et 12 hommes ayant des activités de loisirs différentes.

```r
freqdata <- read.csv('https://marieetienne.github.io/datasets/activite_FC.csv') %>% select(freqC, Activite, Sexe) %>% 
  mutate(Activite = factor(Activite, levels = c('Natation', 'Pilates', 'Pétanque'))) %>% arrange(Activite, Sexe)
freqdata %>% ggplot() +  aes(x= Sexe, y = freqC)+  geom_boxplot(aes(fill = Sexe, col = Sexe), alpha = 0.5) + theme(legend.position = 'none')  + 
  geom_jitter( size=0.8, alpha=0.7, width = 0.15, aes(col = as.factor(Sexe)))
```

--
<p class="question"> Les hommes et les femmes ont-ils la même fréquence cardiaque au repos ?</p>

---
template: intro
## Etude de la fréquence cardiaque au repos

En plus du sexe, il est possible que le type d'activité  influence également la fréquence cardiaque

On souhaite répondre aux questions suivantes:

* Y a t il un effet du sexe  sur la fréquence cardiaque ? 
]

* L'effet de l'activite sur la fréquence cardiaque est-il le même pour les hommes et les femmes ? 
]
---
template: intro
## Cadre général du modèle d'analyse de la variance à 2 facteurs

On étudie le lien entre  
- une variable quantitative notée `$Y$` (la fréquence cardiaque),

Les données peuvent être visualisées à l'aide d'un boxplot.

--
On souhaite répondre aux questions suivantes:

* Y a t il un effet du facteur 1 ?]

* Y a t il un effet du facteur 2 ?

* Y a t il une interaction entre ces deux facteurs ?
]

---
name: model
# Le modèle d'analyse de la variance à 2 facteurs

--
## Graphiquement

Une visualisation graphique du modèle d'analyse de la variance.

Dans l'hypothèse ou le sexe ou l'activité expliquent une partie de la variabilité de la fréquence cardiaque, comment imagine-t-on le processus aléatoire qui a conduit à nos données ?

---
template: model
# Le modèle d'analyse de la variance à 2 facteurs

## Graphiquement

Une visualisation graphique du modèle d'analyse de la variance.

Dans l'hypothèse ou le sexe ou l'activité expliquent une partie de la variabilité de la fréquence cardiaque, comment imagine-t-on le processus aléatoire qui a conduit à nos données ?

---

---
count: false

---
count: false

---
count:false

---

---
count: false

---
count: false

---
template: model

Lequel de ces mécanismes est le plus crédible au vu des donées ?

---
template: model

## Version régulière du modèle Mcomp

`$$\class{alea}{Y_{ijk}} = \class{rouge}{\mu_{ij}} +\class{alea}{E_{ijk}},\quad \class{alea}{E_{ijk}}\overset{ind}{\sim}\mathcal{N}(0, \class{rouge}{\sigma^2}),$$`
avec 
- `$i=1,\ldots,I$` le niveau du facteur 1,
- `$j=1,\ldots,J$` le niveau du facteur 2,
- `$k= 1,\ldots, n_{ij}$` le numéro de l'individu dans le groupe `$(i,j)$`, 
- `$n_{ij}$` le nombre d'individus dans le groupe `$(i,j)$` et `$n=\sum_ j\sum_i n_{ij}$` le nombre total d'individus,
- `$\class{rouge}{\mu_{ij}}$` le comportement moyen du groupe `$(i,j)$`,
- `$\class{rouge}{\sigma^2}$` la variance commune à tous les groupes.

### Nombre de paramètres du modèle

- `$I*J$` paramètres de moyenne  `$(\class{rouge}{\mu_{11}, \mu_{12}, \ldots,  \mu_{1J}, \ldots \mu_{IJ}})$`; 
- 1 paramètre de variance `$\class{rouge}{\sigma^2}$`

---
template: model

## Version régulière du modèle Mcomp sur l'exemple 1
`$$\class{alea}{Y_{ijk}} = \class{rouge}{\mu_{ij}} +\class{alea}{E_{ijk}},\quad \class{alea}{E_{ijk}}\overset{ind}{\sim}\mathcal{N}(0, \class{rouge}{\sigma^2}),$$`
avec 
- `$i=1,\ldots,3$` le niveau du facteur Activite,
- `$j=1,\ldots,2$` le niveau du facteur Sexe,
- `$k= 1,\ldots, 4$` le numéro de l'individu dans le groupe `$(i,j)$`, 
- `$4$` le nombre d'individus dans le groupe `$(i,j)$` et `$n= 24=\sum_ j\sum_i n_{ij}$` le nombre total d'individus,
- `$\class{rouge}{\mu_{ij}}$` le comportement moyen du groupe `$(i,j)$`,
- `$\class{rouge}{\sigma^2}$` la variance commune à tous les groupes.

`$I=3, J= 2$`

### Nombre de paramètres
- 6 paramètres de moyenne
- 1 paramètre de variance

---
template: model
## Version singulière (détaillée) du modèle du modèle Mcomp

`$$\class{alea}{Y_{ijk}} = \class{rouge}{\mu} + \class{rouge}{\alpha_i} + \class{rouge}{\beta_j} + \class{rouge}{\gamma_{ij}} +\class{alea}{E_{ijk}},\quad \class{alea}{E_{ijk}}\overset{ind}{\sim}\mathcal{N}(0, \class{rouge}{\sigma^2}),$$`
avec 
- `$i=1,\ldots,I$` le niveau du facteur 1,
- `$j=1,\ldots,J$` le niveau du facteur 2,
- `$k= 1,\ldots, n_{ij}$` le numéro de l'individu dans le groupe `$(i,j)$`, 
- `$\class{rouge}{\mu}$`  le comportement moyen de référence
- `$\class{rouge}{\alpha_i}$`  l'effet différentiel du niveau `$i$` 
- `$\class{rouge}{\beta_j}$`  l'effet différentiel du niveau `$j$` 
- `$\class{rouge}{\gamma_{ij}}$`  l'effet différentiel de la rencontre des niveaux `$(i,j),$` (terme d'interaction)
- `$\class{rouge}{\sigma^2}$` la variance commune à tous les groupes.

### Nombre de paramètres du modèle

- `$1 + I + J + I*J$` paramètres de moyenne  
- 1 paramètre de variance `$\class{rouge}{\sigma^2}$`

#### La version dans les logiciels et qui permet de séparer les effets des différens facteurs.

---
template: model
## Lien entre les deux versions du même modèle

<table style="width:100%">
  <tr>
    <th>Groupe</th>
    <th>V. régulière</th>
    <th>V. singulière</th>
  </tr>
  <tr>
    <td>1</td>
    <td> `$\mu_{11}$` </td>
    <td> `$\mu +\alpha_1 + \beta_1 + \gamma_{11}$` </td>
  </tr>
  <tr>
    <td>2</td>
    <td> `$\mu_{12}$` </td>
    <td> `$\mu +\alpha_1 + \beta_2 + \gamma_{12}$` </td>
  </tr>
    <tr>
    <td> </td>
    <td>  </td>
    <td>   </td>
  </tr>
  </tr>
    <tr>
    <td> </td>
    <td>  </td>
    <td>   </td>
  </tr>
  </tr>
    <tr>
    <td> </td>
    <td>  </td>
    <td>   </td>
  </tr>
    <tr>
    <td>I</td>
    <td> `$\mu_{IJ}$` </td>
    <td> `$\mu +\alpha_I +\beta_J +\gamma_{IJ}$` </td>
  </tr>
</table>

<a class=care> Problème </a>  Le problème est identique au modèle d'analyse de la variance, le modèle sous cette forme est  <a style="font-weight:400;"> indéterminé</a>.

#### Solution : ajouter des contraintes
Par défaut dans R :
`$\alpha_1=0,$`  `$\beta_1=0,$` et `$\gamma_{1j}=0$` pour tous les `$j$` et `$\gamma_{i1}=0$` pour tous les `$i$`.

Ce qui correspond à `$1 + 1 + J + (I-1)=I  + J +1$` contraintes.

#### Nombre de paramètres 
- On a donc `$1 + I + J + IJ$` paramètres dont `$I  + J +1$` sont contraints, il reste donc
`$IJ$` à estimer pour le comportement moyen
- 1 paramètre de variance `$\sigma^2$`.

---
template: model
## Lien entre les deux versions du même modèle

<table style="width:100%">
  <tr>
    <th>Groupe</th>
    <th>V. régulière</th>
    <th>V. contrainte</th>
  </tr>
  <tr>
    <td>1</td>
    <td> `$\mu_{11}$` </td>
    <td> `$\mu$` </td>
  </tr>
  <tr>
    <td>2</td>
    <td> `$\mu_{12}$` </td>
    <td> `$\mu + \beta_2$` </td>
  </tr>
    <tr>
    <td> </td>
    <td>  </td>
    <td>   </td>
  </tr>
  </tr>
    <tr>
    <td> </td>
    <td>  </td>
    <td>   </td>
  </tr>
  </tr>
    <tr>
    <td> </td>
    <td>  </td>
    <td>   </td>
  </tr>
    <tr>
    <td>I</td>
    <td> `$\mu_{IJ}$` </td>
    <td> `$\mu +\alpha_I +\beta_J +\gamma_{IJ}$` </td>
  </tr>
</table>

---
template: model
## Que signifient ces contraintes (à partir de l'exemple)

* Quelle est la fréquence moyennes des femmes nageuses ?

Les femmes nageuses forment le groupe  `$(1,1)$` : en moyenne leur fréquence cardiaque est donnée par 
`$$\mu +\alpha_1 +\beta_1 + \gamma_{11} = \mu.$$`

* Quelle est la fréquence moyennes des hommes nageurs ?

Les hommes nageurs forment le groupe  `$(1,2)$` : En moyenne leur fréquence cardiaque dans le modèle est 
`$$\mu + \alpha_1 + \beta_2 +\gamma_{12} = \mu + \beta_2.$$`
`$\beta_2$` est donc l'effet différentiel du sexe `$2$` (les hommes)   par rapport au sexe de référence `$1$` (les femmes)   <a class=care> pour la natation</a>.

* Quelle est la fréquence cardiaque moyenne des femmes pratiquant la pétanque ?

Les femmes pratiquant la pétanque  forment le groupe  `$(3,1)$` : En moyenne leur fréquence cardiaque dans le modèle est

`$$\mu + \alpha_3 + \beta_1 +\gamma_{31} = \mu + \alpha_3.$$`

`$\alpha_3$` est donc l'effet différentiel de l'activité 3 (la pétanque) par rapport à l'activité de référence 1 (la natation) pour le sexe 1 (les femmes).

Les valeurs des termes  `$\gamma_{ij}$` sont difficiles à interpréter directement, on en reparlera dans les tests.

---
template: model
## Version contrainte du modèle Mcomp

`$$\class{alea}{Y_{ijk}} = \class{rouge}{\mu} + \class{rouge}{\alpha_i} + \class{rouge}{\beta_j} + \class{rouge}{\gamma_{ij}} +\class{alea}{E_{ijk}},\quad \class{alea}{E_{ijk}}\overset{ind}{\sim}\mathcal{N}(0, \class{rouge}{\sigma^2}),$$`
avec 
- `$i=1,\ldots,I$` le niveau du facteur 1,
- `$j=1,\ldots,J$` le niveau du facteur 2,
- `$k= 1,\ldots, n_{ij}$` le numéro de l'individu dans le groupe `$(i,j)$`, 
- `$\class{rouge}{\mu}$`  le comportement moyen de référence
- `$\class{rouge}{\alpha_i}$`  l'effet différentiel du niveau `$i$`  pour `$i= \class{rouge}{2},\ldots, I$`
- `$\class{rouge}{\beta_j}$`  l'effet différentiel du niveau `$j$` pour `$j = \class{rouge}{2},\ldots, J$`,
- `$\class{rouge}{\gamma_{ij}}$`  l'effet différentiel de la rencontre des niveaux `$(i,j)$`, pour `$\class{rouge}{i\ne1}$` et `$\class{rouge}{j\ne1}$` (terme d'interaction)
- `$\class{rouge}{\sigma^2}$` la variance commune à tous les groupes.

### Nombre de paramètres du modèle

- `$I*J$` paramètres de moyenne  non contraints
- 1 paramètre de variance `$\class{rouge}{\sigma^2}$`

#### La version dans les logiciels et qui permet de séparer les effets des différens facteurs.

---
name: parametre
# Estimation des paramètres
--

###  Estimation des paramètres de moyenne  `$\class{rouge}{(\hat{\mu}, \hat{\alpha_2},\ldots,\hat{\alpha_I}, \hat{\beta_2}, \ldots, \hat{\beta_J}, \hat{\gamma_22}, \ldots,\hat{\gamma_{IJ}})}$`

Par minimisation de la somme des carrés résiduels.

`$$RSS(\mu,\alpha_i,\beta_j,\gamma_{ij}) =  \sum_{i=1}^{I} \sum_{j=1}^{J} \sum_{k=1}^{n_{ij}} \left \lbrace {y_{ijk}} - (  {\mu} + {\alpha_i} + {\beta_j} + {\gamma_{ij}}) \right\rbrace ^2$$`
--

### Estimation de la variance  `$\class{rouge}{\sigma^2}$`

`$$\class{fixe}{\hat{\sigma}}^2 =\frac{1}{n-IJ} \class{fixe}{RSS_{obs}},$$`
`$$\class{fixe}{RSS_{obs}} =  \sum_{i=1}^{I} \sum_{j=1}^{J} \sum_{k=1}^{n_{ij}} \left \lbrace  \class{fixe}{y_k} - (\class{fixe}{\hat{\mu}} + \class{fixe}{\hat{\alpha}_i} + \class{fixe}{\hat{\beta}_j } + \class{fixe}{\hat{\gamma}_{ij}}   ) \right\rbrace ^2$$`

---
name: modcomp
# Test du modèle complet

---
template: modcomp
## Rappel exemple fréquence cardiaque (exemple 1)

On a mesuré la fréquence cardiaque de 12 femmes et 12 hommes pratiquant des activités différentes

--
<p class="question"> Le sexe ou le type d'activité  sont ils liés à la fréquence cardiaque au repos ?</p>

---
template: modcomp
## Sous forme de comparaison de modèle

--
<p class="question"> Le modèle Mcomp est il plus pertinent que le modèle M0 ?</p>

---
template: modcomp
## Hypothèses du test

On va donc opposer une hypothèse de travail `$H_0$` contre une hypothèse alternative `$H_1$`. `$H_0$` peut donc prendre différentes formes:

`$$\begin{align} 
H_0 & =\left \lbrace \mbox{Pas de différence entre les différents groupes }\right\rbrace\\
    & =\left \lbrace  \mbox{pour tout }i, \alpha_i =0, \mbox{pour tout }j,\beta_j=0  \mbox{ et pour tout }(i,j),\gamma_{ij}=0   \right\rbrace\\
    & =\left \lbrace  M_{comp} \mbox{ est équivalent à } M0 \right\rbrace.
\end{align}$$`

`$H_1$` prend les formes équivalentes suivantes

`$$\begin{align} 
H_1 & =\left \lbrace \mbox{Au moins 1 groupe est différent des autres}\right\rbrace\\
    & =\left \lbrace  \mbox{Il existe un }i, \alpha_i  \ne 0  \mbox{ ou un }j, \beta_j \ne 0 \mbox{ ou un } (i,j), \gamma_{ij} \ne 0  \right\rbrace\\
    & =\left \lbrace  M_{comp} \mbox{ est préférable à } M0 \right\rbrace.
\end{align}$$`

Sous `$H_0$`, 
`$$F= \frac{\frac{SS_{M_{comp}}}{IJ-1}}{\frac{RSS}{n-IJ}} \underset{H_0}{\sim}\mathcal{F}(IJ-1, n-IJ)$$`

---
template: modcomp
## Loi de la statistique de test sous `$H_0$` - graphiquement

Sous `$H_0$` la loi de distribution de `$F$` est

---

---
count: false

---
count: false

---
name: variance_decomposition
# Test des effets de chacun des facteurs : décomposition de SSM

---
template: variance_decomposition
## Visualisation graphique de la décomposition de la variance

<br>

- `$RSS_0$` est schématisée par le rectangle  ci dessous.
- La partie rouge correspond à `$RSS$`, la partie blanche à `$RSS_0 -RSS$`, ce que l'on gagne à utilise rle modèle complet

<br>

---
count:false
template: variance_decomposition
## Visualisation graphique de la décomposition de la variance

<br>

- `$RSS_0$` est schématisée par le rectangle  ci dessous.
- La partie rouge correspond à `$RSS$`, la différence en blanc  à `$RSS_0 -RSS$`, ce que l'on gagne à utilise rle modèle complet

<br>

---
count:false
template: variance_decomposition
## Visualisation graphique de la décomposition de la variance

<br>

- `$RSS_0$` est schématisée par le rectangle  ci dessous.

- La partie  beige  `$SSM= RSS_0 -RSS$`, ce que l'on gagne à utiliser le modèle complet plutôt que `$M_0$`.

On souhaite décomposer cette variabilité

---
template: variance_decomposition
## Comment mesurer la variabilité capturée par chaque facteur ?

On veut décomposer la variabilité attribuée au modèle selon les différentes sources.

Considérons les différents modèles possibles

`$$\begin{align}
M_0\ :\ Y_{ijk} &= \mu + E_{ijk}\\
M_{1}\ :\ Y_{ijk} &= \mu + \alpha_i + E_{ijk}\\
M_{2}\ :\ Y_{ijk} &= \mu + \beta_j + E_{ijk}\\
M_{12}\ :\ Y_{ijk} &= \mu + \alpha_i + \beta_j + E_{ijk}\\
M_{comp}\ :\ Y_{ijk} &= \mu + \alpha_i + \beta_j +\gamma_{ij} + E_{ijk}\\
\end{align}$$`

* `$RSS_2 - RSS_{12}$` mesure ce que l'on gagne à utiliser le modèle `$M_{12}$` plutôt que `$M_2$`, c'est l'effet du facteur 1 après s'être ajusté au facteur 1

* `$RSS_1 - RSS_{12}$` mesure ce que l'on gagne à utiliser le modèle `$M_{12}$` plutôt que `$M_1$`, c'est l'effet du facteur 2 après s'être ajusté au facteur 2

* `$RSS_{12} - RSS_{comp}$` mesure ce que l'on gagne à utiliser le modèle `$M_{comp}$` plutôt que `$M_{12}$`, c'est l'effet de l'interaction, les effets principaux étant pris en compte.

---
template: variance_decomposition
## Comment mesurer la variabilité capturée par chaque facteur ?

### Réduction de variabilité

`$$\begin{align}
SSM & = RSS_0 - RSS_{comp}\\
SSM & = \underbrace{RSS_0 - RSS_{1}}_{R(\alpha\vert \mu )} + \underbrace{{RSS_1 - RSS_{12}}}_{R(\beta\vert \mu, \alpha )} + \underbrace{RSS_{12} - RSS_{comp}}_{R(\gamma\vert \mu, \alpha,\beta )} \\
SSM & = \underbrace{RSS_0 - RSS_{2}}_{R(\beta\vert \mu )} + \underbrace{{RSS_2 - RSS_{12}}}_{R(\alpha\vert \mu, \beta )} + \underbrace{RSS_{12} - RSS_{comp}}_{R(\gamma\vert \mu, \alpha,\beta )} \\
\end{align}$$`

On fait un choix pour les parties confondues

---
template: variance_decomposition

## Décomposition de type 2

But : Symétriser le rôle des différents facteurs

---
template: variance_decomposition
## Test de type 2

<table style="width:100%">
  <tr>
    <th>Source</th>
    <th>H0</th>
    <th>SS</th>
    <th>Df</th>
    <th>F</th>
  </tr>
  <tr>
    <th> F1 </th>
    <th> M2 et M12 sont équivalents </th>
    <th> `$R(\alpha\vert \beta, \mu)$` </th>
    <th> `$I-1$` </th>
    <th> `$$\frac{\frac{R(\alpha\vert\beta, \mu)}{I-1}}{\frac{RSS}{n-IJ}}$$` </th>
</tr>
<tr>
    <th> F2 </th>
    <th> M1 et M12 sont équivalents </th>
    <th> `$R(\beta\vert\alpha, \mu)$` </th>
    <th> `$J-1$` </th>
    <th> `$$\frac{\frac{R(\beta\vert \alpha, \mu)}{J-1}}{\frac{RSS}{n-IJ}}$$` </th>
</tr>
<tr>
    <th> Inter </th>
    <th> M12 et Mcomp sont équivalents </th>
    <th> `$R(\gamma\vert \beta, \alpha, \mu)$` </th>
    <th> `$(J-1)(J-1)$` </th>
    <th> `$$\frac{\frac{R(\gamma\vert \beta, \alpha, \mu)}{(I-1)(J-1)}}{\frac{RSS}{n-IJ}}$$` </th>
</tr>
</table>

---
template: variance_decomposition
## L'exemple des fréquences cardiaque

Il s'agit du cas particulier `$n_{ij}$` sont tous égaux. Le plan d'expérience est dit **équilibré**, pas de variabilité confondue.

---
template: variance_decomposition
## L'exemple des fréquences cardiaque

On examine d'abord l'interaction

```r
library(car)
Mcomp <- lm(freqC ~ Sexe + Activite_Fact + Sexe:Activite_Fact, data = freqdata)
Anova(Mcomp, type = 2)
```

```
Anova Table (Type II tests)

Response: freqC
                   Sum Sq Df F value    Pr(>F)    
Sexe                 7.04  1  2.0199   0.17234    
Activite_Fact      410.58  2 58.8884 1.265e-08 ***
Sexe:Activite_Fact  42.58  2  6.1076   0.00945 ** 
Residuals           62.75 18                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

--
 Si l'interaction n'avait pas été significative, on aurait ajuster à nouveau le modèle sans interaction et tester l'effet de chacun des facteurs principaux.

---
name: comp_moyenne
# Comparaison des différents groupes

## Moyennes ajustées

C'est la moyenne qu'on observerait dans un plan d'expérience équilibré

`$$\tilde{\mu}_{i.} = \mu + \alpha_i + \frac{1}{J}\sum_{j=1}^J (\beta_j +\gamma_{ij}).$$`
`$$\tilde{\mu}_{.j} = \mu + \beta_j + \frac{1}{I}\sum_{i=1}^I (\alpha_i +\gamma_{ij}).$$`

`$$\tilde{\mu}_{1.}-\tilde{\mu}_{2.} =  \alpha_1  -\alpha_{2}+ \frac{1}{J}\sum_{j=1}^J (\gamma_{1j} -\gamma_{2j}).$$`

---
template: comp_moyenne

## Moyennes ajustées pour le sexe

```r
library(emmeans)
m_ajust_sexe <- emmeans(Mcomp, "Sexe")
pairs(m_ajust_sexe)
```

```
 contrast estimate    SE df t.ratio p.value
 F - M       -1.08 0.762 18  -1.421  0.1723

Results are averaged over the levels of: Activite_Fact 
```

```r
plot(m_ajust_sexe)
```

---

```r
m_ajust_Activite <- emmeans(Mcomp, "Activite_Fact")
summary(pairs(m_ajust_Activite))
```

```
 contrast            estimate    SE df t.ratio p.value
 Natation - Pilates     -4.75 0.934 18  -5.088  0.0002
 Natation - Pétanque   -10.12 0.934 18 -10.846  <.0001
 Pilates - Pétanque     -5.38 0.934 18  -5.758  0.0001

Results are averaged over the levels of: Sexe 
P value adjustment: tukey method for comparing a family of 3 estimates 
```

---

---

# Ce qu'il faut absolument savoir pour les TDs

- savoir reconnaître une situtation d'analyse de la variance à deux facteurs

--
  * On cherche à comprendre le lien entre une variable d'intérêt `$Y$` et deux variables explicatives qualitatives.
  
--

- savoir écrire le modèle d'analyse de la variance à deux facteurs

--
- Tester l'existence de l'effet d'au moins un des facteurs

--
- Tester l'existence de l'effet d'un facteur en particulier