R : des premiers pas au modèle linéaire

# R : des premiers pas au modèle linéaire
## Modèle linéaire
### Marie-Pierre Etienne
### <a href="https://marieetienne.github.io/premierspas/" class="uri">https://marieetienne.github.io/premierspas/</a>
### Mars 2021

---

---
template: intro

## Des resources pour bien commencer

* <img src =../resources/figs/lm_agro.jpg alt = "Modele lineaire" width = "50" height ="75">    Daudin [Dau15],  [Le modèle linéaire et ses extensions](https://www.editions-ellipses.fr/accueil/3168-statistique-le-modele-lineaire-et-ses-extensions-modele-lineaire-general-modele-lineaire-generalise-modele-mixte-plans-d-experiences-niveau-c-9782340009141.html) (disponible en ligne [ici](https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwjp0vTc67vvAhUNtRoKHQOkD-EQFjABegQIAxAD&url=https%3A%2F%2Fhal.archives-ouvertes.fr%2Fhal-01246652&usg=AOvVaw18GmAe00cx3BqkaIQYTChm) )

* Une [fiche résumé ](https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwjQjv6N7LvvAhWjyIUKHdZCC-UQFjAAegQIBRAD&url=https%3A%2F%2Fwww.math.univ-toulouse.fr%2F~besse%2FWikistat%2Fpdf%2Fst-m-inf-rgls-modlin.pdf&usg=AOvVaw3xIQrAYKSrLMXmyzLNP2cO)

* <img src =../resources/figs/lmr2.png alt = "Modele lineaire" width = "50" height ="75">    Faraway [Far14], [Linear Models with R](https://julianfaraway.github.io/faraway/LMR/)

* La [page du cours](https://marieetienne.github.io/linearmodel.html) que j'ai fait à distance cette année dans un module de Master en écologie

---
template: intro

Le modèle linéaire regroupe un ensemble de méthodes (et de modèles statistiques) qui se traitent de la même manière d'un point de vue mathématique.

* Les modèles d'analyse de la variances
* les modèles de régression simple et multiple
* les modèles d'analyse de la covariance

Dans toutes ses approches, on cherche à caractériser comment les variations d'une grandeur d'intérêt (notée génériquement `$Y$`) s'expliquent par différentes variables explicatives.

---
name: anova1

# Analyse de la variance à 1 facteur

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

On a mesuré la fréquence cardiaque de 20 femmes et 20 hommes.

```r
library(coursesdata) #remotes::install_github('marieetienne/coursesdata')
data(freqdata)
```

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

---
template: anova1
## 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: anova1
## 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=2 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: anova1
## Le modèle d'analyse de la variance à 1 facteur

### Graphiquement

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

---

<img src="modelelineaire_files/figure-html/anova_versiongraphique_proba2_M0-1.png" width="40%" style="display: block; margin: auto;" />
---
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>

---
name:model 
# Analyse de la variance à 1 facteur 
## Le modèle
---
template: model

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

`$$Y_{ik} \overset{ind}{\sim}\mathcal{N}(\mu_i, \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,
- `$\mu_i$` le comportement moyen du groupe `$i$`,
- `$\sigma^2$` la variance commune à tous les groupes.

**Une écriture équivalente **

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

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

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

---
template: model

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

`$$Y_{ik} \overset{ind}{\sim}\mathcal{N}(\mu_i, \sigma^2),$$`
avec  `$I=2$` et  la convention `$i=1$` pour les femmes et `$i=2$` pour les hommes.

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

```
## # A tibble: 2 x 2
##   Sexe      n
##   <chr> <int>
## 1 F        20
## 2 M        20
```
Ainsi 
- `$k=1, \ldots, n_i$`, avec `$n_1=20$` et `$n_2=20$` et `$n=40$`.
- `$\mu_1$` la fréquence cardiaque moyenne des femmes et `$\mu_2$` celle des hommes. 
- `$\sigma^2$` la variance commune à tous les groupes.

**Nombre de paramètres**
- 2 paramètres de moyenne
- 1 paramètre de variance

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

`$$Y_{ik} \overset{ind}{\sim}\mathcal{N}(\mu + \alpha_i, \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$` )
- `$\mu$` un comportement moyen de référence,
- `$\alpha_i$` un effet différentiel du groupe `$i$`par rapport à la référence,
- `$\sigma^2$` la variance commune à tous les groupes.

**Une écriture équivalente **

`$$Y_{ik} = \mu + \alpha_{i} + E_{ik}, \quad E_{ik}\overset{ind}{\sim}\mathcal{N}(0, \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

`$$Y_{ik} \overset{ind}{\sim}\mathcal{N}(\mu + \alpha_i, \sigma^2),$$`
avec 
- `$I=2$` et  la convention `$i=1$` pour les femmes et `$i=2$` pour les hommes.
- `$k=1, \ldots, n_i$`, avec `$n_1=20$` et `$n_2=20$`
- `$\mu$` la fréquence cardiaque moyenne de référence <a class="care"> Référence à définir </a> 
- `$\alpha_1$`, l'effet différentiel d'être une femme par rapport à la différence et `$\alpha_2$` l'effet différentiel d'être un homme.
- `$\sigma^2$` la variance commune à tous les groupes.

**Nombre de paramètres**
- 3 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

`$$Y_{ik} = \mu + \alpha_{i} + E_{ik}, \quad E_{ik}\overset{ind}{\sim}\mathcal{N}(0, \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)
- `$\mu$` un comportement moyen de référence <a class=care> c'est celui du groupe 1</a>
- `$\alpha_1=0$`  la contrainte choisie, `$\alpha_1$` n'est plus un paramètre d'intérêt.
- `$\alpha_i, i\geq 2$` un effet différentiel du groupe `$i$`par rapport à la référence(i.e. le groupe 1)
- `$\sigma^2$` la variance commune à tous les groupes.

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

- `$I$` paramètres de moyenne `$(\mu,  \alpha_2, \ldots, \alpha_I)$`,
- 1 paramètre de variance `$\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

`$$Y_{ik} \overset{ind}{\sim}\mathcal{N}(\mu + \alpha_i, \sigma^2),$$`
avec 
- `$I=2$` et  la convention `$i=1$` pour les femmes et `$i=2$` pour les hommes.
- `$k=1, \ldots, n_i$`, avec `$n_1=20$` et `$n_2=20$`
- `$\mu$` la fréquence cardiaque moyenne des femmesei 
- `$\alpha_1$`, l'effet différentiel d'être une femme par rapport à la différence et `$\alpha_2$` l'effet différentiel d'être un homme.
- `$\sigma^2$` la variance commune à tous les groupes.

**Nombre de paramètres**
- 3 paramètres de moyenne
- 1 paramètre de variance

---
name: parametre
# Analyse de la variance à 1 facteur 
## 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 au choi 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 `$\bf{\theta}$`, `$\mu, \alpha, \dots$`.

On veut en donner une estimation à partir des données observées ( `$\bf{y}$`). Les estimations des paramètres sont notées par les mêmes lettres décorées d'un chapeau `$\bf{\hat{\theta}}, \hat{\mu}, \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 `$\theta$` sera noté `$T$`, l'estimateur de `$\mu$` sera noté `$M$` et l'estimateur de `$\alpha_1$` sera noté `$A_1$`.

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

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

**Estimation**

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

**Estimateur**

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

dont on connaît la loi ([cf cours master pour plus de détails](https://marieetienne.github.io/linear_model/anova.html#38) )

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

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

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

**Estimation**

- `$\mu$` représente le comportement moyen du groupe 1 (groupe de référence) : `$\hat{\mu} = \frac{1}{n_1}\sum_{k=1}^{n_1} y_{1k}$`
--

- `$\mu + \alpha_i = \mu_i$`, donc  `$\alpha_i = \mu_i -\mu$` ,représente l'effet différentiel du groupe i par rapport au groupe de référence.  `$\hat{\alpha_i} = \frac{1}{n_i}\sum_{k=1}^{n_i} y_{ik} - \hat{\mu}$`

**Estimateur **

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

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

L'estimateur de `$\alpha_i$` est donc défini par `$A_i= \frac{1}{n_i}\sum_{k=1}^{n_i} Y_{ik} - \frac{1}{n_1}\sum_{k=1}^{n_1} Y_{1k}.$`

`$$A_i \sim\mathcal{N}\left(\alpha_i, \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 **

`$$\hat{\mu} =\frac{1}{20}\sum_{k=1}^{20} y_{1k}, \quad \hat{\alpha}_2 =\frac{1}{20}\sum_{k=1}^{20} y_{2k} - \frac{1}{20}\sum_{k=1}^{20} y_{1k}$$`

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

```
## # A tibble: 2 x 3
##   Sexe      m     n
##   <chr> <dbl> <int>
## 1 F      84.3    20
## 2 M      84.6    20
```

```r
resM1 <- freqdata %>% group_by(Sexe) %>% summarise(m = mean(freqC), n= n()) %>% 
  mutate(ecart = m -lag(m))
```

`$\hat{\mu} =$` 84.3 et `$\hat{\alpha}_2=$` 0.25.

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

**Ce qu'on doit trouver**

`$\hat{\mu} =$` 84.3 et `$\hat{\alpha}_2=$` 0.25.

**Estimer avec R**

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

```
##             Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)    84.30  0.9025913 93.3977539 1.649393e-46
## SexeM           0.25  1.2764569  0.1958546 8.457675e-01
```

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

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

```
## 
## Call:
## lm(formula = freqC ~ Sexe, data = freqdata)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.550 -2.737 -0.425  2.513  8.700 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  84.3000     0.9026  93.398   <2e-16 ***
## SexeM         0.2500     1.2765   0.196    0.846    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.037 on 38 degrees of freedom
## Multiple R-squared:  0.001008,	Adjusted R-squared:  -0.02528 
## F-statistic: 0.03836 on 1 and 38 DF,  p-value: 0.8458
```

---
name: modcomp
# Analyse de la variance à 1 facteur
## Test du modèle complet

---
template: modcomp

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

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

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

On a mesuré la fréquence cardiaque de 20 femmes et 20 hommes.

--
<p class="question"> Les hommes et les femmes ont-ils la même fréquence cardiaque au repos ?</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 = Variabilité résiduelle dans le modèle `$M_0$`**

`$$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 `$M_1$`**

`$$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 `$M_1$`**

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

<p class=care> Si la variabilité expliquée par le modèle `$M1$` est grande, on gagne à considérer des populations 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 %>%  mutate(m= mean(freqC))
```
]
 
.panel2-rss-user[

```
##    freqC Sport Sexe      m
## 1     82     1    M 84.425
## 2     82     1    F 84.425
## 3     77     1    M 84.425
## 4     79     1    F 84.425
## 5     81     1    M 84.425
## 6     83     1    F 84.425
## 7     76     1    M 84.425
## 8     80     1    F 84.425
## 9     80     2    M 84.425
## 10    83     2    F 84.425
## 11    83     2    M 84.425
## 12    81     2    F 84.425
## 13    81     2    M 84.425
## 14    84     2    F 84.425
## 15    80     2    M 84.425
## 16    82     2    F 84.425
## 17    88     3    M 84.425
## 18    81     3    F 84.425
## 19    87     3    M 84.425
## 20    83     3    F 84.425
## 21    86     3    M 84.425
## 22    83     3    F 84.425
## 23    87     3    M 84.425
## 24    83     3    F 84.425
## 25    85     4    M 84.425
## 26    85     4    F 84.425
## 27    84     4    M 84.425
## 28    86     4    F 84.425
## 29    87     4    M 84.425
## 30    87     4    F 84.425
## 31    85     4    M 84.425
## 32    85     4    F 84.425
## 33    92     5    M 84.425
## 34    88     5    F 84.425
## 35    91     5    M 84.425
## 36    93     5    F 84.425
## 37    90     5    M 84.425
## 38    88     5    F 84.425
## 39    89     5    M 84.425
## 40    90     5    F 84.425
```
]

---
count: false

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

```
##    freqC Sport Sexe      m   Ehat
## 1     82     1    M 84.425 -2.425
## 2     82     1    F 84.425 -2.425
## 3     77     1    M 84.425 -7.425
## 4     79     1    F 84.425 -5.425
## 5     81     1    M 84.425 -3.425
## 6     83     1    F 84.425 -1.425
## 7     76     1    M 84.425 -8.425
## 8     80     1    F 84.425 -4.425
## 9     80     2    M 84.425 -4.425
## 10    83     2    F 84.425 -1.425
## 11    83     2    M 84.425 -1.425
## 12    81     2    F 84.425 -3.425
## 13    81     2    M 84.425 -3.425
## 14    84     2    F 84.425 -0.425
## 15    80     2    M 84.425 -4.425
## 16    82     2    F 84.425 -2.425
## 17    88     3    M 84.425  3.575
## 18    81     3    F 84.425 -3.425
## 19    87     3    M 84.425  2.575
## 20    83     3    F 84.425 -1.425
## 21    86     3    M 84.425  1.575
## 22    83     3    F 84.425 -1.425
## 23    87     3    M 84.425  2.575
## 24    83     3    F 84.425 -1.425
## 25    85     4    M 84.425  0.575
## 26    85     4    F 84.425  0.575
## 27    84     4    M 84.425 -0.425
## 28    86     4    F 84.425  1.575
## 29    87     4    M 84.425  2.575
## 30    87     4    F 84.425  2.575
## 31    85     4    M 84.425  0.575
## 32    85     4    F 84.425  0.575
## 33    92     5    M 84.425  7.575
## 34    88     5    F 84.425  3.575
## 35    91     5    M 84.425  6.575
## 36    93     5    F 84.425  8.575
## 37    90     5    M 84.425  5.575
## 38    88     5    F 84.425  3.575
## 39    89     5    M 84.425  4.575
## 40    90     5    F 84.425  5.575
```
]

---
count: false

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

```
##    freqC Sport Sexe      m   Ehat     E2hat
## 1     82     1    M 84.425 -2.425  5.880625
## 2     82     1    F 84.425 -2.425  5.880625
## 3     77     1    M 84.425 -7.425 55.130625
## 4     79     1    F 84.425 -5.425 29.430625
## 5     81     1    M 84.425 -3.425 11.730625
## 6     83     1    F 84.425 -1.425  2.030625
## 7     76     1    M 84.425 -8.425 70.980625
## 8     80     1    F 84.425 -4.425 19.580625
## 9     80     2    M 84.425 -4.425 19.580625
## 10    83     2    F 84.425 -1.425  2.030625
## 11    83     2    M 84.425 -1.425  2.030625
## 12    81     2    F 84.425 -3.425 11.730625
## 13    81     2    M 84.425 -3.425 11.730625
## 14    84     2    F 84.425 -0.425  0.180625
## 15    80     2    M 84.425 -4.425 19.580625
## 16    82     2    F 84.425 -2.425  5.880625
## 17    88     3    M 84.425  3.575 12.780625
## 18    81     3    F 84.425 -3.425 11.730625
## 19    87     3    M 84.425  2.575  6.630625
## 20    83     3    F 84.425 -1.425  2.030625
## 21    86     3    M 84.425  1.575  2.480625
## 22    83     3    F 84.425 -1.425  2.030625
## 23    87     3    M 84.425  2.575  6.630625
## 24    83     3    F 84.425 -1.425  2.030625
## 25    85     4    M 84.425  0.575  0.330625
## 26    85     4    F 84.425  0.575  0.330625
## 27    84     4    M 84.425 -0.425  0.180625
## 28    86     4    F 84.425  1.575  2.480625
## 29    87     4    M 84.425  2.575  6.630625
## 30    87     4    F 84.425  2.575  6.630625
## 31    85     4    M 84.425  0.575  0.330625
## 32    85     4    F 84.425  0.575  0.330625
## 33    92     5    M 84.425  7.575 57.380625
## 34    88     5    F 84.425  3.575 12.780625
## 35    91     5    M 84.425  6.575 43.230625
## 36    93     5    F 84.425  8.575 73.530625
## 37    90     5    M 84.425  5.575 31.080625
## 38    88     5    F 84.425  3.575 12.780625
## 39    89     5    M 84.425  4.575 20.930625
## 40    90     5    F 84.425  5.575 31.080625
```
]

---
count: false

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

]

---
count: false

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

```
## # A tibble: 40 x 4
## # Groups:   Sexe [2]
##    freqC Sport Sexe      m
##    <dbl> <int> <chr> <dbl>
##  1    82     1 M      84.6
##  2    82     1 F      84.3
##  3    77     1 M      84.6
##  4    79     1 F      84.3
##  5    81     1 M      84.6
##  6    83     1 F      84.3
##  7    76     1 M      84.6
##  8    80     1 F      84.3
##  9    80     2 M      84.6
## 10    83     2 F      84.3
## # … with 30 more rows
```
]

---
count: false

```
## # A tibble: 40 x 5
## # Groups:   Sexe [2]
##    freqC Sport Sexe      m  Ehat
##    <dbl> <int> <chr> <dbl> <dbl>
##  1    82     1 M      84.6 -2.55
##  2    82     1 F      84.3 -2.30
##  3    77     1 M      84.6 -7.55
##  4    79     1 F      84.3 -5.30
##  5    81     1 M      84.6 -3.55
##  6    83     1 F      84.3 -1.30
##  7    76     1 M      84.6 -8.55
##  8    80     1 F      84.3 -4.30
##  9    80     2 M      84.6 -4.55
## 10    83     2 F      84.3 -1.30
## # … with 30 more rows
```
]

---
count: false

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

```
## # A tibble: 40 x 6
## # Groups:   Sexe [2]
##    freqC Sport Sexe      m  Ehat E2hat
##    <dbl> <int> <chr> <dbl> <dbl> <dbl>
##  1    82     1 M      84.6 -2.55  6.50
##  2    82     1 F      84.3 -2.30  5.29
##  3    77     1 M      84.6 -7.55 57.0 
##  4    79     1 F      84.3 -5.30 28.1 
##  5    81     1 M      84.6 -3.55 12.6 
##  6    83     1 F      84.3 -1.30  1.69
##  7    76     1 M      84.6 -8.55 73.1 
##  8    80     1 F      84.3 -4.30 18.5 
##  9    80     2 M      84.6 -4.55 20.7 
## 10    83     2 F      84.3 -1.30  1.69
## # … with 30 more rows
```
]

---
count: false

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

```
## # A tibble: 40 x 6
##    freqC Sport Sexe      m  Ehat E2hat
##    <dbl> <int> <chr> <dbl> <dbl> <dbl>
##  1    82     1 M      84.6 -2.55  6.50
##  2    82     1 F      84.3 -2.30  5.29
##  3    77     1 M      84.6 -7.55 57.0 
##  4    79     1 F      84.3 -5.30 28.1 
##  5    81     1 M      84.6 -3.55 12.6 
##  6    83     1 F      84.3 -1.30  1.69
##  7    76     1 M      84.6 -8.55 73.1 
##  8    80     1 F      84.3 -4.30 18.5 
##  9    80     2 M      84.6 -4.55 20.7 
## 10    83     2 F      84.3 -1.30  1.69
## # … with 30 more rows
```
]

---
count: false

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

]

---
count: false

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

```
## [1] 619.775
```

```
## [1] 619.15
```
]

---
count: false

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

```
## [1] 619.775
```

```
## [1] 619.15
```

```
## [1] 0.625
```
]

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

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

```
## Analysis of Variance Table
## 
## Model 1: freqC ~ 1
## Model 2: freqC ~ Sexe
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     39 619.77                           
## 2     38 619.15  1     0.625 0.0384 0.8458
```

---
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$`,

`$$\frac{SS_{M1}}{\sigma^2} \underset{H_0}{\sim}\chi^2(I-1)$$`

--
Mais `$\sigma^2$` est inconnu, on a envie de le remplacer par son estimateur `$RSS / (n-I).$`

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

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

### Proabilité critique

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

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

Grossièrement dit, 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)
```

<img src="modelelineaire_files/figure-html/pvalue_graphique-ex_FC-1.png" width="60%" style="display: block; margin: auto;" />
---
template: modcomp
### Exercice 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>

```r
install.packages('remotes')
remotes::install_github('marieetienne/coursesdata')
data(bats)
```

---
name: test_param
# Analyse de la variance à 1 facteur
## Test sur un paramètre

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

`$$Y_{ik} = \mu +\alpha_i + E_{ik}, \quad i=1,\ldots I,\ k=1,\ldots, n_i.$$`
Par défaut, R propose pour chaque coefficient le test de l'hypothèse `$H_0\ :\ \left\lbrace\mbox{ le coefficient est nul } \right\rbrace$`.

Que signifie ce test ?

- pour `$\mu$`
- pour `$\alpha_i$`

<a class=care> Attention, le sens de ce test dépend  de la contrainte choisie </p>

Pour plus de détail sur les tests sur des paramètres se reporter au [cours en ligne](https://marieetienne.github.io/linearmodel.html)
* extension à plus que 1 facteur

---
# Analyse de la variance à 1 facteur
## 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

---
# Analyse de la variance à 1 facteur
## 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 Sexe).

### 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),

Sur un  paramètre  pvaut une valeur donnée (par défaut R teste systématiquement si un paramètre vaut 0)
- on doit vérifier que les hypothèses de modélisation ne sont pas trop fausses.

---

---
template: reg
## Etude de la pollution au SO2

On a mesuré pour 41 villes américaines, la pollution au SO2 ainsi que la population dans la ville

```r
library(coursesdata)
data(usdata)
```

--
<p class="question"> La taille d'une ville est-elle liée à la pollution en SO2 ?</p>

---
template: reg
## Cadre général du modèle de régression simple

On étudie le lien entre  
- une variable quantitative notée `$Y$` (l'indicateur de SO2),
- et une variable quantitative `$x$`.

Les données peuvent être visualisées à l'aide d'un nuage de points.

--
<p class="question"> La variable x permet elle d'expliquer la variabilité de la variable Y ?</p>

---
name: model
# Regression linéaire 
## Le modèle de régression simple

---
template: model

### Graphiquement

Une visualisation graphique du modèle d'analyse de régression simple

Comment imagine-t-on le processus aléatoire qui a conduit à nos données ?

---

---
count: false

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

`$$Y_{k} = \beta_0 +\beta_1 x_{k}  +E_{k},\quad E_{k}\overset{ind}{\sim}\mathcal{N}(0, \sigma^2),$$`
avec 
- `$x_k$` la valeur de la variable explicative pour l'observation `$k$`, 
- `$k=1,\ldots,n$` le numéro d'individu, `$n$` le nombre total d'individus,
- `$\beta_0$` l'ordonnée à l'origine, 
- `$\beta_1$` la pente de la droite, mesure de l'effet de la variable `$x$`

- `$\sigma^2$` la variance.

### Une écriture équivalente

`$$Y_{k} \overset{ind}{\sim}\mathcal{N}(\beta_0+\beta_1 x_k, \sigma^2).$$`

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

- `$2$` paramètres de moyenne  `$(\beta_0, \beta_1)$`; 
- 1 paramètre de variance `$\sigma^2$`

---
template: model

## Le modèle de régression simple sur l'exemple de la pollution.

`$$Y_{k} = \beta_0 +\beta_1 x_{k}  +E_{k},\quad E_{k}\overset{ind}{\sim}\mathcal{N}(0, \sigma^2),$$`
avec 
- `$x_k$` la population dans la ville `$k$`, 
- `$k=1,\ldots,n$` le numéro de la ville, `$n=41$`,
- `$\beta_0$` l'ordonnée à l'origine, 
- `$\beta_1$` la pente de la droite, mesure de l'effet de la population sur la pollution.

- `$\sigma^2$` la variance.

### Nombre de paramètres du modèle
- 2 paramètres de moyennes
- 1 paramètre de variance

---
name: parametre
# Regression linéaire
## Estimation des paramètres

---
template: parametre

### Version mathématique

hors programme !

---
template: parametre
### Version pratique

```r
Mpop <- lm(SO2~pop, data = usdata)
summary(Mpop)
```

```
## 
## Call:
## lm(formula = SO2 ~ pop, data = usdata)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.545 -14.456  -4.019  11.019  72.549 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 17.868316   4.713844   3.791 0.000509 ***
## pop          0.020014   0.005644   3.546 0.001035 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 20.67 on 39 degrees of freedom
## Multiple R-squared:  0.2438,	Adjusted R-squared:  0.2244 
## F-statistic: 12.57 on 1 and 39 DF,  p-value: 0.001035
```

---
name: prediction
# Regeression linéaire
## Prediction avec un modèle de régression simple

---
template:  prediction
Il est fréquent d'utiliser un modèle de régression pour prédire.

### Prédiction de la valeur moyenne pour un `$x$` particulier

**Valeur moyenne attendue pour `$y$` pour un `$x$` donné**
`$\beta_0+\beta_1 x.$`

**Valeur moyenne prédite pour `$y$` pour un `$x$` donné**
`$\hat{\beta}_0+\hat{\beta}_1 x.$`

```r
predict(Mpop, newdata=data.frame(pop=333))
```

```
##        1 
## 24.53284
```

### Intervalle de confiance pour la valeur moyenne prédite  pour `$y$` pour un `$x$` donné
`$\hat{\beta}_0+\hat{\beta}_1 x.$`

```r
predict(Mpop, newdata=data.frame(pop=333), interval = 'confidence')
```

```
##        fit      lwr      upr
## 1 24.53284 17.28453 31.78115
```

---
template:  prediction
Il est fréquent d'utiliser un modèle de régression pour prédire.

### Prédiction de la valeur possible `$y$` pour un  `$x$` particulier

```r
predict(Mpop, newdata=data.frame(pop=333), interval = 'prediction')
```

```
##        fit       lwr      upr
## 1 24.53284 -17.90225 66.96793
```

---
template:  prediction
### Sur l'exemple de la polution

**Intervalle de confiance pour le comportement moyen**

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

---
count: false
template:  prediction
### Sur l'exemple de la polution

**Intervalle de confiance pour le comportement moyen**

```
## Warning in predict.lm(Mpop, interval = "prediction", level = 0.95): predictions on current data refer to _future_ responses
```

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

---
name: model
# Modèle de régression multiple

---
template: model

## Le modèle de régression multiple

Plusieurs variables sont potentiellement liées à la pollution en SO2

- temp : Average temperature in Fahrenheit
- manuf : No. of companies employing more than 20 employees
- pop : Population in thousands
- wind : Average annual wind speed in miles/hour
- precip : annual precipitation height in inches
- days : No. of days of precipitation

<a class=question> Quelles sont les variables liées à la pollution en SO2 ? </a>

---
template: model

## Le modèle de régression multiple

`$$Y_{k} = \beta_0 +\beta_1 x_{k}^{1}  + \beta_2 x_{k}^{2} + \ldots +  \beta_p x_{k}^{p}  +  E_{k},\quad E_{k}\overset{ind}{\sim}\mathcal{N}(0, \sigma^2),$$`
avec 
- `$x_{k}^{l}$` la valeur de la variable explicative `$l$` pour l'observation `$k$`, 
- `$k=1,\ldots,n$` le numéro d'individu, `$n$` le nombre total d'individus,
- `$\beta_0$` l'ordonnée à l'origine, 
- `$\beta_l$` l'effet de la variable `$X^{l}$` sur la variable à expliquer,
- `$\sigma^2$` la variance.

### Une écriture équivalente

`$$Y_{k} \overset{ind}{\sim}\mathcal{N}(\beta_0 +\beta_1 x_{k}^{1}  + \beta_2 x_{k}^{2} + \ldots +  \beta_p x_{k}^{p} , \sigma^2).$$`

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

- `$l+1$` paramètres de moyenne  `$(\beta_0, \beta_1, \ldots, \beta_l)$`; 
- 1 paramètre de variance `$\sigma^2$`

---
template: model

### Sur l'exemple de la pollution

```r
Mcomp <- lm(SO2 ~ temp + manuf + pop + wind + precip + days, data = usdata) 
#Mcomp <- lm(SO2 ~ . - City , data = usdata) # toutes les variables sauf City
model.matrix(Mcomp) %>% head(n = 3)
```

```
##   (Intercept) temp manuf pop wind precip days
## 1           1 70.3   213 582  6.0   7.05   36
## 2           1 61.0    91 132  8.2  48.52  100
## 3           1 56.7   453 716  8.7  20.66   67
```
---
name: parametre
# Regression linéaire multiple
## Estimation des paramètres

---
template: parametre
### Estimation des paramètres du modèle version mathématique

hors programme

### Estimation des paramètres du modèle version pratique

```r
Mcomp<- lm(SO2~ temp +  manuf + pop + wind +  precip + days, data = usdata)
summary(Mcomp)
```

```
## 
## Call:
## lm(formula = SO2 ~ temp + manuf + pop + wind + precip + days, 
##     data = usdata)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.004  -8.542  -0.991   5.758  48.758 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 111.72848   47.31810   2.361 0.024087 *  
## temp         -1.26794    0.62118  -2.041 0.049056 *  
## manuf         0.06492    0.01575   4.122 0.000228 ***
## pop          -0.03928    0.01513  -2.595 0.013846 *  
## wind         -3.18137    1.81502  -1.753 0.088650 .  
## precip        0.51236    0.36276   1.412 0.166918    
## days         -0.05205    0.16201  -0.321 0.749972    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.64 on 34 degrees of freedom
## Multiple R-squared:  0.6695,	Adjusted R-squared:  0.6112 
## F-statistic: 11.48 on 6 and 34 DF,  p-value: 5.419e-07
```

---
template: parametre
### Estimation des paramètres du modèle version mathématique

hors programme

### Estimation des paramètres du modèle version mathématique
On ne peut prendre en compte que certaines  variables

Exemple

```r
M12<- lm(SO2~pop + wind, data = usdata)
summary(M12)
```

```
## 
## Call:
## lm(formula = SO2 ~ pop + wind, data = usdata)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.363 -14.562  -4.336  11.019  72.794 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 19.486939  22.161866   0.879  0.38477   
## pop          0.020107   0.005851   3.436  0.00144 **
## wind        -0.177390   2.371746  -0.075  0.94077   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 20.94 on 38 degrees of freedom
## Multiple R-squared:  0.2439,	Adjusted R-squared:  0.2041 
## F-statistic:  6.13 on 2 and 38 DF,  p-value: 0.004928
```

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

---
template: modcomp
## Pollution

<p class="question"> La pollution en SO2 dans les villes américaines est elles liées à l'une au moins des variables caractérisiques des villes ?</p>

On va à la pêche ....

---
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{Auncune variable n'est liée à la pollution en SO2}\right\rbrace\\
    & =\left \lbrace  \mbox{pour tout }p\geq 1, \beta_p =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 variable est liée à la pollution en SO2}\right\rbrace\\
    & =\left \lbrace  \mbox{Il existe un } p, \beta_p \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}}}{l}}{\frac{RSS}{n-l-1}} \underset{H_0}{\sim}\mathcal{F}(l, n-l-1)$$`

---
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: test_variable
# Regression multiple 
## Test de l'effet des variables

---
template: test_variable

### Test sur les paramètres

Tester la nullité du paramètre `$\beta_l$` revient à tester si la variable `$x^{l}$` et la variable `$Y$` sont liées.

Ce test est similaire  au test de comparaison entre le modèle complet et le modèle complet privé de la variable `$x^{l}$`.

---
count: false
template: test_variable

### Equivalence des tests sur l'exemple de la pollution

```
## Loading required package: carData
```

```
## 
## Attaching package: 'car'
```

```
## The following object is masked from 'package:dplyr':
## 
##     recode
```

```
## The following object is masked from 'package:purrr':
## 
##     some
```

```r
summary(Mcomp)$coefficients
```

```
##                 Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept) 111.72848064 47.31810073  2.361221 0.0240867374
## temp         -1.26794109  0.62117952 -2.041183 0.0490557189
## manuf         0.06491817  0.01574825  4.122245 0.0002277862
## pop          -0.03927674  0.01513274 -2.595482 0.0138461970
## wind         -3.18136579  1.81501910 -1.752800 0.0886503978
## precip        0.51235896  0.36275507  1.412410 0.1669175999
## days         -0.05205019  0.16201386 -0.321270 0.7499724652
```

```r
Mcomp_l <- lm(SO2 ~  temp + manuf +  wind + precip + days, data = usdata) 
anova(Mcomp_l, Mcomp)
```

```
## Analysis of Variance Table
## 
## Model 1: SO2 ~ temp + manuf + wind + precip + days
## Model 2: SO2 ~ temp + manuf + pop + wind + precip + days
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1     35 8726.3                              
## 2     34 7283.3  1    1443.1 6.7365 0.01385 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

### Lien entre les statistiques  de tests

```r
res <- summary(Mcomp)$coefficients
res[,"t value"]^2
```

```
## (Intercept)        temp       manuf         pop        wind      precip 
##   5.5753634   4.1664282  16.9929075   6.7365251   3.0723084   1.9949026 
##        days 
##   0.1032144
```

---
template: test_variable
## Vigilance sur l'interprétation des tests

```r
summary(Mcomp)$coefficients
```

```r
GGally::ggpairs(usdata, columns = 2:8)
```

---

# References

Cornillon, P., A. Guyader, F. Husson, et al. (2018). _R pour la
statistique et la science des données_. Presses universitaires de
Rennes.

Daudin, J. (2015). _Le modèle linéaire et ses extensions-Modèle
linéaire général, modèle linéaire généralisé, modèle mixte, plans
d'expériences (Niveau C)_. Edition Ellipses.

Faraway, J. J. (2014). _Linear models with R_. CRC press.

Hutcheon, J. M., J. A. W. Kirsch, and T. Garland (2002). "A Comparative
Analysis of Brain Size in Relation to Foraging Ecology and Phylogeny in
the Chiroptera". In: _Brain, Behavior and Evolution_ 60.3, pp. 165-180.
ISSN: 1421-9743. DOI:
[10.1159/000065938](https://doi.org/10.1159%2F000065938). URL:
[http://dx.doi.org/10.1159/000065938](http://dx.doi.org/10.1159/000065938).

R Core Team (2017). _R: A Language and Environment for Statistical
Computing_. R Foundation for Statistical Computing. Vienna, Austria.
URL: [https://www.R-project.org/](https://www.R-project.org/).

Wickham, H. and G. Grolemund (2016). _R for data science: import, tidy,
transform, visualize, and model data_. " O'Reilly Media, Inc.".