Analyse de la variance

# Analyse de la variance
## à 1 facteur
### Marie-Pierre Etienne
### <a href="https://github.com/marieetienne" class="uri">https://github.com/marieetienne</a>
### 2020/09/11 (updated: 2023-09-08)

---

---
template: intro
## Un exemple jouet à but pédagogique

On a mesuré la fréquence cardiaque de 12 femmes et 12 hommes qui pratique des activités sportives différentes (Natation, Pilate, Pétanque).

```r
freqdata <- read.table('https://marieetienne.github.io/datasets/activite_FC.csv', sep = ",", header = TRUE )
freqdata %>% ggplot() +  aes(x= Activite, y = freqC) +  geom_boxplot( alpha = 0.5) + theme(legend.position = 'none')  + 
  geom_jitter(  alpha=0.7, width = 0.15)
```

--
<p class="question"> L'activité physique est-elle liée à  fréquence cardiaque au repos ?</p>

---
template: intro
## Cas d'étude Chauve souris

- Des espèces avec différents régimes alimentaires.
- L'écho localisation (capacité auditive) est un sens essentiel pour les insectivores.

Hutcheon, Kirsch, and Garland [HKG02] ont recensé pour 63 espèces de chauves souris l'espèce, le régime alimentaire, le clade (groupe dans l'arbre phylogénétique), différentes variables morphologiques dont le volume du cerveau dédié à l'audition (AUD)

--
<p class="question"> Le volume de la partie auditive du cerveau est il lié au régime alimentaire ?</p>

---
template: intro
## Cadre général du modèle d'analyse de la variance à 1 facteur
On étudie le lien entre  
- variable quantitative notée `$Y$` (la fréquence cardiaque, volume de cerveau),
- et une variable qualitative (un facteur) pourvant prendre `$I$` modalités (I=3 pour l'exemple 1 et I=4 dans l'exemple 2).

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

--
<p class="question"> La variable d'intérêt Y  est-elle en moyenne différente selon les différentes modalités ?</p>

---
template: intro

## Graphiquement

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

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
count: false

<p class="question"> Au vu de nos données, le modèle M1 est-il plus pertinent que le modèle M0 ?</p>

---
# Plan du cours

## Formaliser les modèles M0 et M1

## Estimer les paramètres inconnus

## Proposer une méthode pour décider entre M0 et M1

## Si M1, comparer les comportements moyens des groupes 2 à 2

---
name: model
# Formaliser les modèles M0 et M1

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

`$$\class{alea}{Y_{ik}} \overset{ind}{\sim}\mathcal{N}(\class{rouge}{\mu_i}, \class{rouge}{\sigma^2}),$$`
avec 
- `$i=1,\ldots,I$` le numéro du groupe,
- `$k= 1,\ldots, n_i$` le numéro de l'individu dans le groupe `$i$`, 
- `$n_i$` le nombre d'individus dans le groupe `$i$` et `$n=\sum_i n_i$` le nombre total d'individus,
- `$\class{rouge}{\mu_i}$` le comportement moyen du groupe `$i$`,
- `$\class{rouge}{\sigma^2}$` la variance commune à tous les groupes.

### Une écriture équivalente

`$$\class{alea}{Y_{ik}} = \class{rouge}{\mu_{i}} + \class{alea}{E_{ik}}, \quad \class{alea}{E_{ik}}\overset{ind}{\sim}\mathcal{N}(0, \class{rouge}{\sigma^2}).$$`

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

- `$I$` paramètres de moyenne  `$(\class{rouge}{\mu_1}, \class{rouge}{\mu_2}, \ldots, \class{rouge}{\mu_I})$`; 
- 1 paramètre de variance `$\class{rouge}{\sigma^2}$`

---
template: model

## Version régulière du modèle M1 sur l'exemple 1

`$$\class{alea}{Y_{ik}} \overset{ind}{\sim}\mathcal{N}(\class{rouge}{\mu_i}, \class{rouge}{\sigma^2}),$$`
avec  `$I=3$` et  la convention `$i=1$` pour la natation et `$i=2$` pour le pilate et `$i=3$` pour la pétanque.

```r
freqdata %>%group_by(Activite) %>% summarise(n= n())
```

```
# A tibble: 3 × 2
  Activite     n
  <chr>    <int>
1 Natation     8
2 Pétanque     8
3 Pilates      8
```

- `$k=1, \ldots, n_i$`, avec `$n_1=8$` et `$n_2=8$` et `$n_3=8$` et `$n = 24$`.
- `$\class{rouge}{\mu_1}$` la fréquence cardiaque moyenne des nageurs/nageuses  et `$\class{rouge}{\mu_2}$` des pratiquant.e.s de pilates et `$\class{rouge}{\mu_3}$` celle des joueurs/joueuses de pétanque. 
- `$\class{rouge}{\sigma^2}$` la variance commune à tous les groupes.

### Nombre de paramètres
- 3 paramètres de moyenne et  1 paramètre de variance

---
template: model
## Version singulière du modèle du modèle M1

`$$\class{alea}{Y_{ik}} \overset{ind}{\sim}\mathcal{N}(\class{rouge}{\mu} + \class{rouge}{\alpha_i}, \class{rouge}{\sigma^2}),$$`
avec 
- `$i=1,\ldots,I$` le numéro du groupe,
- `$k= 1,\ldots, n_i$` le numéro de l'individu dans le groupe `$i$`, ( `$\sum_i n_i=n$` )
- `$\class{rouge}{\mu}$` un comportement moyen de référence,
- `$\class{rouge}{\alpha_i}$` un effet différentiel du groupe `$i$`par rapport à la référence,
- `$\class{rouge}{\sigma^2}$` la variance commune à tous les groupes.

### Une écriture équivalente

`$$\class{alea}{Y_{ik}} = \class{rouge}{\mu} + \class{rouge}{\alpha_{i}} + \class{alea}{E_{ik}}, \quad \class{alea}{E_{ik}}\overset{ind}{\sim}\mathcal{N}(0, \class{rouge}{\sigma^2}).$$`

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

- `$I+1$` paramètres de moyenne `$(\mu, \alpha_1, \alpha_2, \ldots, \alpha_I)$`,
- 1 paramètre de variance `$\sigma^2.$`

#### La version dans les logiciels et celle qui se généralise à plusieurs facteurs.

---
template: model
## Version singulière du modèle du modèle M1 sur l'exemple 1

`$$\class{alea}{Y_{ik}} \overset{ind}{\sim}\mathcal{N}(\class{rouge}{\mu} + \class{rouge}{\alpha_i}, \class{rouge}{\sigma^2}),$$`
avec 
- `$I=3$` et  la convention `$i=1$` pour la natation et `$i=2$` pour le pilate et `$i=3$` pour la pétanque.
- `$k=1, \ldots, n_i$`, avec `$n_1=8,$`  `$n_2=8$` et `$n_3=8$`, 
- `$\class{rouge}{\mu}$` la fréquence cardiaque moyenne de référence <a class="care"> Référence à définir </a> 
- `$\class{rouge}{\alpha_1}$`, l'effet différentiel de la pratique de la natation  par rapport à la référence, `$\class{rouge}{\alpha_2}$` l'effet différentiel de la pratique du pilate et `$\class{rouge}{\alpha_3}$` l'effet différentiel de la pratique de la pétanque  par rapport à la référence.
- `$\class{rouge}{\sigma^2}$` la variance commune à tous les groupes.

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

---
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_1$` </td>
    <td> `$\mu +\alpha_1$` </td>
  </tr>
  <tr>
    <td>2</td>
    <td> `$\mu_2$` </td>
    <td> `$\mu +\alpha_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_I$` </td>
    <td> `$\mu +\alpha_I$` </td>
  </tr>
</table>

<a class=care> Problème </a> : Version régulière mal définie. Le modèle dans cette version est dit <a style="font-weight:400;"> indéterminé</a>.

--
#### Exemple
Si `$\mu_1 =10,\ \mu_2=12,$` dans la forme singulière peut correspondre à

- `$\mu =10,\ \alpha_1=0, \ \alpha_2=2,$`
- ou `$\mu =0,\ \alpha_1=10, \ \alpha_2=12,$`
- ou `$\mu =11,\ \alpha_1=-1, \ \alpha_2=1,$`
- ou `$\mu =15,\ \alpha_1=-5, \ \alpha_2=-3,$`
- `$\ldots$`

--
<a class=care> Solution </a> : Choisir une contrainte. Par défaut dans R, `$\alpha_1=0$`. La référence est le comportement du groupe 1.

---
template: model
## Version singulière du modèle du modèle M1 en intégrant la contrainte par défaut

`$$\class{alea}{Y_{ik}} = \class{rouge}{\mu} + \class{rouge}{\alpha_{i}} + \class{alea}{E_{ik}}, \quad \class{alea}{E_{ik}}\overset{ind}{\sim}\mathcal{N}(0, \class{rouge}{\sigma^2}).$$`

avec 
- `$i=1,\ldots,I$` le numéro du groupe,
- `$k= 1,\ldots, n_i$` le numéro de l'individu dans le groupe `$i$`, ($\sum_i n_i = n$ nombre total d'individus)
- `$\class{rouge}{\mu}$` un comportement moyen de référence <a class=care> c'est celui du groupe 1</a>
- `$\class{rouge}{\alpha_1}=0$`  la contrainte choisie, `$\alpha_1$` n'est plus un paramètre d'intérêt.
- `$\class{rouge}{\alpha_i}, i\geq 2$` un effet différentiel du groupe `$i$`par rapport à la référence(i.e. le groupe 1)
- `$\class{rouge}{\sigma^2}$` la variance commune à tous les groupes.

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

- `$I$` paramètres de moyenne `$(\class{rouge}{\mu},  \class{rouge}{\alpha_2}, \ldots, \class{rouge}{\alpha_I})$`,
- 1 paramètre de variance `$\class{rouge}{\sigma^2}.$`

---
template: model
## Version singulière du modèle du modèle M1 sur l'exemple 1  en intégrant la contrainte par défaut

`$$\class{alea}{Y_{ik}} \overset{ind}{\sim}\mathcal{N}(\class{rouge}{\mu} + \class{rouge}{\alpha_i}, \class{rouge}{\sigma^2}),$$`
avec 
- `$I=3$` et  la convention `$i=1$` pour la natation et `$i=2$` pour le pilate et `$i=3$` pour la pétanque.
- `$k=1, \ldots, n_i$`, avec `$n_1=8,$`  `$n_2=8$` et `$n_3=8$`, 
- `$\class{rouge}{\mu}$` la fréquence cardiaque moyenne de référence, c'est-à-dire, avec la contrainte par défaut de R la fréquence cardiaque des nageurs/nageuses. 
- `$\class{rouge}{\alpha_1}$`, l'effet différentiel de la pratique de la natation  par rapport à la référence (donc valant 0), `$\class{rouge}{\alpha_2}$` l'effet différentiel de la pratique du pilate et `$\class{rouge}{\alpha_3}$` l'effet différentiel de la pratique de la pétanque  par rapport à la référence.
- `$\class{rouge}{\sigma^2}$` la variance commune à tous les groupes.

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

---
name: parametre
# Etude des paramètres du modèle

---
template: parametre
## Quelques précisions de notations et de vocabulaires

- Les données observées (les valeurs prises par la variable `$Y$` ) : `$\bf{y}=(y_{11}, y_{12}, \ldots, y_{I n_I}).$`

- En refaisant la même expérience, on obtiendrait d'autres valeurs (liées àaux choix d'individus différents par exemple, ou mesuré un autre jour, ...), on va noter `$Y_{ik}$` la variable aléatoire qui décrit la loi des valeurs possibles pour l'individu k du groupe i.

<br>
<a class=care> les minuscules notent des valeurs, les majuscules des variables aléatoires. </a>
<br>

Les paramètres du modèle  sont inconnus. Ils sont notés par des lettres grecques `$\class{rouge}{\bf{\theta}}$`, `$\class{rouge}{\mu}, \class{rouge}{\alpha}, \dots$`.

On veut en donner une estimation à partir des données observées ( `$\class{fixe}{\bf{y}}$`). Les estimations des paramètres sont notées par les mêmes lettres décorées d'un chapeau `$\class{fixe}{\bf{\hat{\theta}}}, \class{fixe}{\hat{\mu}}, \class{fixe}{\hat{\alpha}}, \dots$`.

L'estimateur d'un paramètre est la variable aléatoire correspondante. Il sera noté en majuscule avec une lettre latine. Par exemple, l'etimateur de `$\class{rouge}{\theta}$` sera noté `$\class{alea}{T}$`, l'estimateur de `$\class{rouge}{\mu}$` sera noté `$\class{alea}{M}$` et l'estimateur de `$\class{rouge}{\alpha_1}$` sera noté `$\class{alea}{A_1}$`.

---
template: parametre
## Estimation des paramètres du modèle version régulière

`$$\class{alea}{Y_{ik}} = \class{rouge}{\mu_i} + \class{alea}{E_{ik}}, \quad \class{alea}{E_{ik}}\overset{ind}{\sim}\mathcal{N}(0,\class{rouge}{\sigma^2}).$$`

### Estimation

`$\class{rouge}{\mu_i}$` représente le comportement moyen de l'ensemble des individus du groupe `$i$`.  On peut l'estimer par `$\class{fixe}{\hat{\mu}_i} = \frac{1}{n_i}\sum_{k=1}^{n_i} \class{fixe}{y_{ik}}= \class{fixe}{y_{i\bullet}}$`

### Estimateur

L'estimateur de `$\class{rouge}{\mu_i}$` est donc défini par `$\class{alea}{M_i}= \frac{1}{n_i}\sum_{k=1}^{n_i} \class{alea}{Y_{ik}},$`  c'est une variable aléatoire

- de loi normale `$\mathcal{N}$`, 
--

- et espérance `$\mathbb{E}(\class{alea}{M_i}) = \class{rouge}{\mu_i}.$` 
--
<a class=care> C'est un estimateur sans biais</a>
--

- de variance `$\mathbb{V}(\class{alea}{M_i}) =  \frac{\class{rouge}{\sigma^2}}{n_i}$`
--
<a class=care> la variance de l'estimateur diminue quand le nombre d'individus dans le groupe augmente.</a>
--

`$$\class{alea}{M_i} \sim\mathcal{N}\left(\class{rouge}{\mu_i}, \frac{\class{rouge}{\sigma^2}}{n_i}\right).$$`

---
template: parametre
## Estimation des paramètres du modèle version singulière avec la contrainte par défaut

`$$\class{alea}{Y_{ik}} = \class{rouge}{\mu} + \class{rouge}{\alpha_i} + \class{alea}{E_{ik}}, \quad \class{alea}{E_{ik}}\overset{ind}{\sim}\mathcal{N}(0,\class{rouge}{\sigma^2}).$$`

### Estimation

- `$\class{rouge}{\mu}$` représente le comportement moyen du groupe 1 (groupe de référence) : `$\class{fixe}{\hat{\mu}} =  \class{fixe}{y_{1\bullet}}$`
--

- `$\class{rouge}{\mu} + \class{rouge}{\alpha_i} = \class{rouge}{\mu_i}$`, donc  `$\class{rouge}{\alpha_i} = \class{rouge}{\mu_i} -\class{rouge}{\mu}$` ,représente l'effet différentiel du groupe i par rapport au groupe de référence.  `$\class{fixe}{\hat{\alpha_i}} = \class{fixe}{y_{i\bullet}} -\class{fixe}{y_{1\bullet}}$`

### Estimateur

L'estimateur de `$\class{alea}{\mu}$` est donc défini par `$\class{alea}{M}= \class{alea}{Y_{1\bullet}} = \frac{1}{n_1}\sum_{k=1}^{n_1} \class{alea}{Y_{1k}}.$`

`$$\class{alea}{M} \sim\mathcal{N}\left(\class{rouge}{\mu}, \frac{\class{rouge}{\sigma^2}}{n_1}\right).$$`

L'estimateur de `$\alpha_i$` est donc défini par `$\class{alea}{A_i}=  \class{alea}{Y_{i\bullet}} -\class{alea}{Y_{1\bullet}}.$`

`$$\class{alea}{A_i} \sim\mathcal{N}\left(\class{rouge}{\alpha_i}, \class{rouge}{\sigma^2}\left(\frac{1}{n_i} + \frac{1}{n_1}\right) \right).$$`

---
template: parametre
## Estimation des paramètres pour l'exemple 1.

### Ce qu'on doit trouver

`$$\class{fixe}{\hat{\mu}} = \class{fixe}{y_{1\bullet}}=\frac{1}{8}\sum_{k=1}^{8} y_{1k}= \frac{y_{11} + y_{12} + \ldots + y_{18}}{8}$$`
`$$\quad \class{fixe}{\hat{\alpha}_2} =\class{fixe}{y_{2\bullet}}-\class{fixe}{y_{1\bullet}}=\frac{1}{8}\sum_{k=1}^{8} \class{fixe}{y_{2k}} - \frac{1}{8}\sum_{k=1}^{8} \class{fixe}{y_{1k}}, \quad \class{fixe}{\hat{\alpha}_3} =\class{fixe}{y_{3\bullet}}-\class{fixe}{y_{1\bullet}}=\frac{1}{8}\sum_{k=1}^{8} \class{fixe}{y_{3k}} - \frac{1}{8}\sum_{k=1}^{8} \class{fixe}{y_{1k}}$$`

```r
freqdata %>% group_by(Activite) %>% summarise(m = mean(freqC), n= n())
```

```
# A tibble: 3 × 3
  Activite     m     n
  <chr>    <dbl> <int>
1 Natation  70       8
2 Pétanque  80.1     8
3 Pilates   74.8     8
```

`$\class{fixe}{\hat{\mu}} =$` 70 , `$\class{fixe}{\hat{\alpha}_2}=$` 10.125 et `$\class{fixe}{\hat{\alpha}_3}=$` 4.75.

---
count: false
template: parametre
## Estimation des paramètres pour l'exemple 1.

### Ce qu'on doit trouver

`$\class{fixe}{\hat{\mu}} =$` 70 et `$\class{fixe}{\hat{\alpha}_2}=$` 10.125et `$\class{fixe}{\hat{\alpha}_3}=$` 4.75.

### Estimer avec R

```r
M1 <- lm(freqC ~ Activite, data = freqdata)
summary(M1)$coefficients
```

```
                 Estimate Std. Error   t value     Pr(>|t|)
(Intercept)        70.000  0.8178623 85.588976 3.355090e-28
ActivitePétanque   10.125  1.1566320  8.753865 1.887487e-08
ActivitePilates     4.750  1.1566320  4.106751 5.033922e-04
```

---
count: false
template: parametre
## Estimation des paramètres pour l'exemple 1.

```r
M1 <- lm(freqC ~ Activite, data = freqdata)
summary(M1)
```

```

Call:
lm(formula = freqC ~ Activite, data = freqdata)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0000 -1.7500 -0.0625  2.0000  3.2500

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       70.0000     0.8179  85.589  < 2e-16 ***
ActivitePétanque  10.1250     1.1566   8.754 1.89e-08 ***
ActivitePilates    4.7500     1.1566   4.107 0.000503 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.313 on 21 degrees of freedom
Multiple R-squared:  0.7851,	Adjusted R-squared:  0.7647 
F-statistic: 38.36 on 2 and 21 DF,  p-value: 9.73e-08
```

---
template: parametre
## Retour sur la méthode d'estimation - Somme des carrés résiduel (RSS)

Naturellement, on a choisi d'estimer  `$\class{rouge}{\mu_i}$` par `$\class{fixe}{y_{i\bullet}}.$`

En fait, `$\class{fixe}{\hat{\mu}_i}$` est choisi pour minimiser `$\class{fixe}{RSS_{M_1, obs}}$`, la somme observée des erreurs au carré  donnée par

`$$\class{fixe}{RSS_{M_1, obs}} = \sum_{i=1}^I \sum_{k=1}^{n_i} (\class{fixe}{y_{ik}} - \class{fixe}{\hat{\mu}_i})^2.$$` 
--

Le modèle d'analyse de la variance 
`$$\class{alea}{Y_{ik}} =\class{rouge}{\mu_i} +\class{alea}{E_{ik}}.$$`
Pour un nouvel individu `$k_0$` dans le groupe `$i$`, le modèle prédit que sa fréquence cardiaque sera `$\class{fixe}{\hat{y}_{ik_0}}=\class{fixe}{\hat{\mu}_i}$`.

`$\class{fixe}{RSS_{M_1,obs}}$` se réécrit dans un cadre général comme

`$$\class{fixe}{RSS_{M_1, obs}} = \sum_{i=1}^I \sum_{k=1}^{n_i} (\class{fixe}{y_{ik}} - \class{fixe}{\hat{y}_{ik}})^2.$$`

`$$\class{alea}{RSS_{M_1}} = \sum_{i=1}^I \sum_{k=1}^{n_i} (\class{alea}{Y_{ik}} - \class{alea}{\hat{Y}_{ik}})^2.$$`

---
template: parametre

## Loi de RSS

`$$\frac{\class{alea}{RSS}}{\class{rouge}{\sigma^2}} \sim \chi^2(n-I).$$`

--
## Estimateur de la variance

$$\class{alea}{S^2} =\frac{1}{n-I} \class{alea}{RSS}, $$
est un <a class=care> estimateur sans bias de  `$\class{rouge}{\sigma^2}$` </a> .

## Estimation de `$\sigma^2$`

`$$\class{fixe}{\hat{\sigma}}^2 =\frac{1}{n-I} \class{fixe}{RSS_{obs}}.$$`

---

```r
*freqdata %>%  group_by(Activite)
```
]
 
.panel2-var_est-user[

```
# A tibble: 24 × 3
# Groups:   Activite [3]
  freqC Sexe  Activite
  <int> <chr> <chr>   
1    72 M     Natation
2    72 F     Natation
3    67 M     Natation
4    69 F     Natation
# … with 20 more rows
```
]

---
count: false

```r
freqdata %>%  group_by(Activite) %>%
* mutate(yhat =  mean(freqC))
```
]
 
.panel2-var_est-user[

```
# A tibble: 24 × 4
# Groups:   Activite [3]
  freqC Sexe  Activite  yhat
  <int> <chr> <chr>    <dbl>
1    72 M     Natation    70
2    72 F     Natation    70
3    67 M     Natation    70
4    69 F     Natation    70
# … with 20 more rows
```
]

---
count: false

```r
freqdata %>%  group_by(Activite) %>%
  mutate(yhat =  mean(freqC)) %>%
* mutate( Ehat = freqC-yhat,Ehat_2 = Ehat^2)
```
]
 
.panel2-var_est-user[

```
# A tibble: 24 × 6
# Groups:   Activite [3]
  freqC Sexe  Activite  yhat  Ehat Ehat_2
  <int> <chr> <chr>    <dbl> <dbl>  <dbl>
1    72 M     Natation    70     2      4
2    72 F     Natation    70     2      4
3    67 M     Natation    70    -3      9
4    69 F     Natation    70    -1      1
# … with 20 more rows
```
]

---
count: false

```r
freqdata %>%  group_by(Activite) %>%
  mutate(yhat =  mean(freqC)) %>%
  mutate( Ehat = freqC-yhat,Ehat_2 = Ehat^2) %>%
* ungroup %>% summarise(RSSobs=sum(Ehat_2)) %>% as.numeric() -> RSSobs
*RSSobs
```
]
 
.panel2-var_est-user[

```
[1] 112.375
```
]

---
count: false

```r
freqdata %>%  group_by(Activite) %>%
  mutate(yhat =  mean(freqC)) %>%
  mutate( Ehat = freqC-yhat,Ehat_2 = Ehat^2) %>%
  ungroup %>% summarise(RSSobs=sum(Ehat_2)) %>% as.numeric() -> RSSobs
RSSobs
*sigma2_hat <- RSSobs / (nrow(freqdata) - 2)
*sigma2_hat
```
]
 
.panel2-var_est-user[

```
[1] 112.375
```

```
[1] 5.107955
```
]

---
count: false

```r
freqdata %>%  group_by(Activite) %>%
  mutate(yhat =  mean(freqC)) %>%
  mutate( Ehat = freqC-yhat,Ehat_2 = Ehat^2) %>%
  ungroup %>% summarise(RSSobs=sum(Ehat_2)) %>% as.numeric() -> RSSobs
RSSobs
sigma2_hat <- RSSobs / (nrow(freqdata) - 2)
sigma2_hat
*sigma_hat <- sqrt(sigma2_hat)
*sigma_hat
```
]
 
.panel2-var_est-user[

```
[1] 112.375
```

```
[1] 5.107955
```

```
[1] 2.260078
```
]

---

### On doit trouver

`$\hat{\sigma}=$` 2.313264.

### Avec R

```r
M1 <- lm(freqC ~ Activite, data = freqdata)
```

```r
summary(M1)$sigma
```

```
[1] 2.313264
```

---
count: false
template: parametre
## Estimation du paramètre de variance sur l'exemple 1

### On doit trouver

`$\hat{\sigma}=$` 2.313264.

```r
summary(M1)
```

```

Call:
lm(formula = freqC ~ Activite, data = freqdata)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0000 -1.7500 -0.0625  2.0000  3.2500

Residual standard error: 2.313 on 21 degrees of freedom
Multiple R-squared:  0.7851,	Adjusted R-squared:  0.7647 
F-statistic: 38.36 on 2 and 21 DF,  p-value: 9.73e-08
```

---
template: parametre
## Visualisation graphique des paramètres estimés

---

---
count: false

---
count: false

---
template: pause

### Méfiez vous

- vérifier les degrés de liberté et le nom des paramètres
- Si quelque chose cloche, vérifier la nature des variables
- On peut transformer une variable en facteur à l'aide de la fonction `as.factor`.

---
name: modcomp
# Décider entre M0 et M1

**Test du modèle complet**

## Rappel exemple fréquence cardiaque (exemple 1)

On a mesuré la fréquence cardiaque de 8 nageurs/nageuses, 8 pratiquant.e.s de pîlate et 8 pratiquant.e.s de pétanque

--
<p class="question"> Y a -t-il un lien entre fréquence cardiaque et activité sportive ?</p>

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

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

Autrement dit

---
template: modcomp

---
template: modcomp
## Test de l'existence d'un effet

<p class="question"> Le comportement moyen de la variable `$Y$` est il différent selon les différentes modalités de la variable explicative ?</p>

Deux modèles possibles
`$$M1 : Y_{ik} = \mu  + \alpha_i + E_{ik},\quad  M0: Y_{ik} = \mu  + E_{ik}$$`

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

---
template: modcomp
## Mesurer la variabilité expliquée par le modèle M1

### Variabilité totale dans les données

`$$RSS_{0, obs} = \sum_{i=1}^{I} \sum_{k=1}^{n_i} (y_{ik} - \bar{y})^2 = \sum_{i=1}^{I} \sum_{k=1}^{n_i} (y_{ik} - \hat{y}_{ik}^{(0)})^2.$$` 
--

### Variabilité résiduelle dans le modèle M1

`$$RSS_{1, obs} = \sum_{i=1}^{I} \sum_{k=1}^{n_i}  ( y_{ik} - (\hat{\mu} +\hat{\alpha}_i))^2 = \sum_{i=1}^{I} \sum_{k=1}^{n_i}  ( y_{ik} - \hat{y}_{ik}^{(1)})^2.$$` 
--

### Variabilité expliquée par le modèle M1

`$$SCM = SS_{M1,obs} =  RSS_{0, obs} - RSS_{1, obs}.$$` 
C'est la quantité de variabilité expliquée par les différence entre les I populations.

<p class=care> Si la variabilité expliquée par le modèle `$M1$` est grande, on gagne à considérer des ppopulations différentes <p>

On a mis en évidence des différences entre les populations.

---
template: modcomp
## Calcul de RSS sur l'exemple 1

A la main :

---

```r
*freqdata %>% group_by(Activite) %>%
* mutate(m= mean(freqC))
```
]
 
.panel2-rss-user[

```
# A tibble: 24 × 4
# Groups:   Activite [3]
  freqC Sexe  Activite     m
  <int> <chr> <chr>    <dbl>
1    72 M     Natation    70
2    72 F     Natation    70
3    67 M     Natation    70
4    69 F     Natation    70
# … with 20 more rows
```
]

---
count: false

```r
freqdata %>% group_by(Activite) %>%
  mutate(m= mean(freqC)) %>%
* mutate( Ehat = freqC - m)
```
]
 
.panel2-rss-user[

```
# A tibble: 24 × 5
# Groups:   Activite [3]
  freqC Sexe  Activite     m  Ehat
  <int> <chr> <chr>    <dbl> <dbl>
1    72 M     Natation    70     2
2    72 F     Natation    70     2
3    67 M     Natation    70    -3
4    69 F     Natation    70    -1
# … with 20 more rows
```
]

---
count: false

```r
freqdata %>% group_by(Activite) %>%
  mutate(m= mean(freqC)) %>%
  mutate( Ehat = freqC - m) %>%
* mutate( E2hat = Ehat^2)
```
]
 
.panel2-rss-user[

```
# A tibble: 24 × 6
# Groups:   Activite [3]
  freqC Sexe  Activite     m  Ehat E2hat
  <int> <chr> <chr>    <dbl> <dbl> <dbl>
1    72 M     Natation    70     2     4
2    72 F     Natation    70     2     4
3    67 M     Natation    70    -3     9
4    69 F     Natation    70    -1     1
# … with 20 more rows
```
]

---
count: false

```r
freqdata %>% group_by(Activite) %>%
  mutate(m= mean(freqC)) %>%
  mutate( Ehat = freqC - m) %>%
  mutate( E2hat = Ehat^2) %>%
* ungroup()
```
]
 
.panel2-rss-user[

```
# A tibble: 24 × 6
  freqC Sexe  Activite     m  Ehat E2hat
  <int> <chr> <chr>    <dbl> <dbl> <dbl>
1    72 M     Natation    70     2     4
2    72 F     Natation    70     2     4
3    67 M     Natation    70    -3     9
4    69 F     Natation    70    -1     1
# … with 20 more rows
```
]

---
count: false

```r
freqdata %>% group_by(Activite) %>%
  mutate(m= mean(freqC)) %>%
  mutate( Ehat = freqC - m) %>%
  mutate( E2hat = Ehat^2) %>%
  ungroup() %>%
* summarise(RSSM1 =sum(E2hat)) %>% as.numeric()
```
]
 
.panel2-rss-user[

```
[1] 112.375
```
]

---
count: false

]

---
count: false

```
   freqC Sexe Activite        m
1     72    M Natation 74.95833
2     72    F Natation 74.95833
3     67    M Natation 74.95833
4     69    F Natation 74.95833
5     71    M Natation 74.95833
6     73    F Natation 74.95833
7     66    M Natation 74.95833
8     70    F Natation 74.95833
9     78    M  Pilates 74.95833
10    71    F  Pilates 74.95833
11    77    M  Pilates 74.95833
12    73    F  Pilates 74.95833
13    76    M  Pilates 74.95833
14    73    F  Pilates 74.95833
15    77    M  Pilates 74.95833
16    73    F  Pilates 74.95833
17    82    M Pétanque 74.95833
18    78    F Pétanque 74.95833
19    81    M Pétanque 74.95833
20    83    F Pétanque 74.95833
21    80    M Pétanque 74.95833
22    78    F Pétanque 74.95833
23    79    M Pétanque 74.95833
24    80    F Pétanque 74.95833
```
]

---
count: false

```r
freqdata %>% group_by(Activite) %>%
  mutate(m= mean(freqC)) %>%
  mutate( Ehat = freqC - m) %>%
  mutate( E2hat = Ehat^2) %>%
  ungroup() %>%
  summarise(RSSM1 =sum(E2hat)) %>% as.numeric()->
  RSSM1
freqdata %>%  mutate(m= mean(freqC)) %>%
* mutate( Ehat = freqC - m)
```
]
 
.panel2-rss-user[

```
   freqC Sexe Activite        m      Ehat
1     72    M Natation 74.95833 -2.958333
2     72    F Natation 74.95833 -2.958333
3     67    M Natation 74.95833 -7.958333
4     69    F Natation 74.95833 -5.958333
5     71    M Natation 74.95833 -3.958333
6     73    F Natation 74.95833 -1.958333
7     66    M Natation 74.95833 -8.958333
8     70    F Natation 74.95833 -4.958333
9     78    M  Pilates 74.95833  3.041667
10    71    F  Pilates 74.95833 -3.958333
11    77    M  Pilates 74.95833  2.041667
12    73    F  Pilates 74.95833 -1.958333
13    76    M  Pilates 74.95833  1.041667
14    73    F  Pilates 74.95833 -1.958333
15    77    M  Pilates 74.95833  2.041667
16    73    F  Pilates 74.95833 -1.958333
17    82    M Pétanque 74.95833  7.041667
18    78    F Pétanque 74.95833  3.041667
19    81    M Pétanque 74.95833  6.041667
20    83    F Pétanque 74.95833  8.041667
21    80    M Pétanque 74.95833  5.041667
22    78    F Pétanque 74.95833  3.041667
23    79    M Pétanque 74.95833  4.041667
24    80    F Pétanque 74.95833  5.041667
```
]

---
count: false

```r
freqdata %>% group_by(Activite) %>%
  mutate(m= mean(freqC)) %>%
  mutate( Ehat = freqC - m) %>%
  mutate( E2hat = Ehat^2) %>%
  ungroup() %>%
  summarise(RSSM1 =sum(E2hat)) %>% as.numeric()->
  RSSM1
freqdata %>%  mutate(m= mean(freqC)) %>%
  mutate( Ehat = freqC - m) %>%
* mutate( E2hat = Ehat^2)
```
]
 
.panel2-rss-user[

```
   freqC Sexe Activite        m      Ehat     E2hat
1     72    M Natation 74.95833 -2.958333  8.751736
2     72    F Natation 74.95833 -2.958333  8.751736
3     67    M Natation 74.95833 -7.958333 63.335069
4     69    F Natation 74.95833 -5.958333 35.501736
5     71    M Natation 74.95833 -3.958333 15.668403
6     73    F Natation 74.95833 -1.958333  3.835069
7     66    M Natation 74.95833 -8.958333 80.251736
8     70    F Natation 74.95833 -4.958333 24.585069
9     78    M  Pilates 74.95833  3.041667  9.251736
10    71    F  Pilates 74.95833 -3.958333 15.668403
11    77    M  Pilates 74.95833  2.041667  4.168403
12    73    F  Pilates 74.95833 -1.958333  3.835069
13    76    M  Pilates 74.95833  1.041667  1.085069
14    73    F  Pilates 74.95833 -1.958333  3.835069
15    77    M  Pilates 74.95833  2.041667  4.168403
16    73    F  Pilates 74.95833 -1.958333  3.835069
17    82    M Pétanque 74.95833  7.041667 49.585069
18    78    F Pétanque 74.95833  3.041667  9.251736
19    81    M Pétanque 74.95833  6.041667 36.501736
20    83    F Pétanque 74.95833  8.041667 64.668403
21    80    M Pétanque 74.95833  5.041667 25.418403
22    78    F Pétanque 74.95833  3.041667  9.251736
23    79    M Pétanque 74.95833  4.041667 16.335069
24    80    F Pétanque 74.95833  5.041667 25.418403
```
]

---
count: false

```r
freqdata %>% group_by(Activite) %>%
  mutate(m= mean(freqC)) %>%
  mutate( Ehat = freqC - m) %>%
  mutate( E2hat = Ehat^2) %>%
  ungroup() %>%
  summarise(RSSM1 =sum(E2hat)) %>% as.numeric()->
  RSSM1
freqdata %>%  mutate(m= mean(freqC)) %>%
  mutate( Ehat = freqC - m) %>%
  mutate( E2hat = Ehat^2) %>%
* summarise(RSSM1 =sum(E2hat)) %>% as.numeric()
```
]
 
.panel2-rss-user[

```
[1] 522.9583
```
]

---
count: false

```r
freqdata %>% group_by(Activite) %>%
  mutate(m= mean(freqC)) %>%
  mutate( Ehat = freqC - m) %>%
  mutate( E2hat = Ehat^2) %>%
  ungroup() %>%
  summarise(RSSM1 =sum(E2hat)) %>% as.numeric()->
  RSSM1
freqdata %>%  mutate(m= mean(freqC)) %>%
  mutate( Ehat = freqC - m) %>%
  mutate( E2hat = Ehat^2) %>%
  summarise(RSSM1 =sum(E2hat)) %>% as.numeric()->
* RSSM0
*SSM <- RSSM0 - RSSM1
```
]
 
.panel2-rss-user[

]

---
count: false

```r
freqdata %>% group_by(Activite) %>%
  mutate(m= mean(freqC)) %>%
  mutate( Ehat = freqC - m) %>%
  mutate( E2hat = Ehat^2) %>%
  ungroup() %>%
  summarise(RSSM1 =sum(E2hat)) %>% as.numeric()->
  RSSM1
freqdata %>%  mutate(m= mean(freqC)) %>%
  mutate( Ehat = freqC - m) %>%
  mutate( E2hat = Ehat^2) %>%
  summarise(RSSM1 =sum(E2hat)) %>% as.numeric()->
  RSSM0
SSM <- RSSM0 - RSSM1
*SSM
```
]
 
.panel2-rss-user[

```
[1] 410.5833
```
]

```
[1] 410.5833
```

---
template: modcomp
## Calcul de RSS sur l'exemple 1

```r
M1 <- lm(freqC ~ Activite, data= freqdata)
M0 <- lm(freqC ~ 1, data= freqdata)
anova(M0, M1)
```

```
Analysis of Variance Table

Model 1: freqC ~ 1
Model 2: freqC ~ Activite
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)    
1     23 522.96                                 
2     21 112.37  2    410.58 38.364 9.73e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---
template: modcomp
## Quand décider si `$SS_{M1, obs}$` est grande ?

--
C'est le rôle du test statistique.

Idée de démarche:

- Décrire quelle gamme de valeurs pourrait prendre SSM dans le cas où il n'y pas de différence entre les groupes 
- Comparer avec ce qu'on constate sur nos données
- 
    - Si c'est comparable, pas de raison de penser qu'il y a des différences entre les groupes,
    - Si c'est surprenant, les groupes sont vraisemblablement différentes.

On veut décrire comment varie `$SS_{M1, obs}$` dans la situation où il n"y a pas de différence entre les groupes.

On travaille donc sous l'hypothèse 
`$$H_0 =\left \lbrace \mbox{Pas de différence entre les groupes }\right\rbrace$$`

`$$H_0 =\left \lbrace  \mbox{pour tout }i, \alpha_i =0  \right\rbrace$$`

`$$H_0 =\left \lbrace  M1 \mbox{ est équivalent à } M0 \right\rbrace.$$`

---
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 groupes }\right\rbrace\\
    & =\left \lbrace  \mbox{pour tout }i, \alpha_i =0  \right\rbrace\\
    & =\left \lbrace  M1 \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  \right\rbrace\\
    & =\left \lbrace  M1 \mbox{ est préférable à } M0 \right\rbrace.
\end{align}$$`

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

Il faut être capable de décrire quelles sont les valeurs possibles de `$SS_{M1, obs}$`, pour ceci il faut connaitre la loi de la variable aléatoire `$SS_{M1}$`.

On peut montrer que sous `$H_0$`,

<p class=question> 
Sous `$H_0$`, 
`$$F= \frac{\frac{SS_{M1}}{I-1}}{\frac{RSS}{n-I}} \underset{H_0}{\sim}\mathcal{F}(I-1, n-I)$$`  
</p>

Intuitivement le numérateur de `$F$` mesure la force d'un signal, le dénominateur un niveau de bruit. Plus la valeur de `$F$` est grande, plus il est raisonnable de penser que le signal capturé par le modèle est fort et donc révélateur d'un véritable effet de la variable explicative.

---
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

---
template: modcomp

## Mesurer le caractère surprenant : Proabilité critique

--
La probabilité critique ( p-value en anglais) est définie par

$$ \mathbb{P}(F> F_{obs} \vert H_0)$$

Intuitivement c'est la probabilité d'observer Fobs ou des valeurs encore plus extrèmes sous `$H_0$`.

On va rejeter `$H_0$` lorsque la probabilité critique est faible (typiquement inférieure à 5%). La valeur de la statistique de test observée est peu compatible avec l'hypothèse `$H_0$`. On ne croit pas à `$H_0$`.

---
template: modcomp
## Déclinaison sur l'exemple Fréquence cardiaque

```r
anova(M0,M1)
```

```
Analysis of Variance Table

<img src="anova_files/figure-html/pvalue_graphique-ex_FC-1.png" width="60%" height="30%" style="display: block; margin: auto;" />
---
template: modcomp
## Sur l'exemple chauve souris

### Pouvez vous répondre à la question

<a class=question> Y a-t-il un effet du régime alimentaire sur le volume de la partie auditive du cerveau ? </a>

---
name: test_param
# Comparer les modalités

---
template: test_param
## Réflexion sur le sens de ce test

On souhaite comparer les différentes modalités : et répondre aux questions

* La fréquence cardiaque des nageurs est-elle différente de la fréquence cardiaque des pratiquants de yoga ?
* La fréquence cardiaque des nageurs est-elle différente de la fréquence cardiaque des joueurs de pétanque ? 
* La fréquence cardiaque des pratiquants de yoga est-elle différente de la fréquence cardiaque des joueurs de pétanque ?

Comment faire ?

Tester

$$H_0 : \left\lbrace (\mu+\alpha_i) = (\mu + \alpha_j) \right\rbrace = \left\lbrace\alpha_i =  \alpha_j \right\rbrace $$
contre 
$$H_1 : \left \lbrace(\mu+\alpha_i \ne \mu + \alpha_j) \right\rbrace = \left\lbrace\alpha_i \ne  \alpha_j \right\rbrace $$
---
template: test_param
## Mise en place du test de comparaison des modalités

Sous `$H_0$`,

`$$T = \frac{Y_{i,\bullet} -Y_{j,\bullet}}{ \hat{\sigma} \sqrt{\frac{1}{n_i} + \frac{1}{n_j}}} \sim \mathcal{T}(DDL_{res})$$`

où `$\mathcal{T}(DDL_{res})$` désigne la loi de Student à `$DDL_{res}$` degrés de liberté. 
--

### Définition de la p-value

---
count: false

---
count: false

---
count: false

---
count: false

---
count: false

---
count: false

---
count: false

---
count: false

---
count: false

---
template: test_param

## Déclinaison sur les exemples

## Comparaison des pratiques sportives

---

```
 Activite emmean    SE df lower.CL upper.CL
 Natation   70.0 0.818 21     68.3     71.7
 Pétanque   80.1 0.818 21     78.4     81.8
 Pilates    74.8 0.818 21     73.0     76.5

Confidence level used: 0.95 
```

```
 contrast            estimate   SE df t.ratio p.value
 Natation - Pétanque   -10.12 1.16 21  -8.754  <.0001
 Natation - Pilates     -4.75 1.16 21  -4.107  0.0015
 Pétanque - Pilates      5.38 1.16 21   4.647  0.0004

P value adjustment: bonferroni method for 3 tests 
```

---

- `$\hat{\sigma}=$` 2.313.
- `$\sqrt{c_j}\hat{\sigma}=$` 0.818,
- `$t_{j,obs}=$` 85.589,
- p-value= 0,

```r
summary(M1)
```

```

Call:
lm(formula = freqC ~ Activite, data = freqdata)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0000 -1.7500 -0.0625  2.0000  3.2500

Residual standard error: 2.313 on 21 degrees of freedom
Multiple R-squared:  0.7851,	Adjusted R-squared:  0.7647 
F-statistic: 38.36 on 2 and 21 DF,  p-value: 9.73e-08
```

---
# Démarche complète d'analyse à partir de l'exemple des chauves souris

Une question initiale :  Le volume auditif dépend il de régime alimentaire. 
- Représenter les données en fonction des questions qu'on se pose.
- Ecrire le modèle 
- Traduire la question en 1 ou plusieurs tests statistiques
- Apporter une réponse concrète

---
# Bilan sur le modèle d'analyse de la variance à 1 facteur

### But 
Etudier le lien entre une variable quantitative (la fréquence cardiaque) et un facteur (le Activite).

## Le modèle 
`$$Y_{ik} = \mu + \alpha_i + E_{ik}$$`
avec `$\mu$` le comportement de référence et `$\alpha_i$` l'effet différentiel du groupe `$i$`  par rapport à la référence.

La référence est définie par la contrainte choisie. R par défaut choisit la contrainte `$\alpha_1=0$`, ce qui place le groupe 1 comme groupe de référence.

## Tests

* Pour comparer des modèles ( pour tester un potentiel effet du facteur),

* Comparer les modalités (trouver quelles sont les modalités différentes les unes des autres)

## Conclusion concrète

---