Régression linéaire

# Régression linéaire
## multiple
### 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 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) #remotes::install_github('MarieEtienne/coursesdata')
data(usdata)
```

---
name: reg_multiple
# Le modèle de régression multiple

## Présentation du problème
.pull-left[
Dans l'exemple de la pollution en SO2 des villes, 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 tibble: 41 × 8
  City            SO2  temp manuf   pop  wind precip  days
  <chr>         <int> <dbl> <int> <int> <dbl>  <dbl> <int>
1 Phoenix          10  70.3   213   582   6     7.05    36
2 Little rock      13  61      91   132   8.2  48.5    100
3 San Francisco    12  56.7   453   716   8.7  20.7     67
4 Denver           17  51.9   454   515   9    13.0     86
5 Hartford         56  49.1   412   158   9    43.4    127
# … with 36 more rows
```

]

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

---
template: reg_multiple

## Le modèle de régression multiple

`$$\class{alea}{Y_{k}} = \class{rouge}{\beta_0} +\class{rouge}{\beta_1} x_{k,1}  + \class{rouge}{\beta_2} x_{k,2} + \ldots +  \class{rouge}{\beta_p} x_{k,p}  +  \class{alea}{E_{k}},\quad \class{alea}{E_{k}}\overset{ind}{\sim}\mathcal{N}(0, \class{rouge}{\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,
- `$\class{rouge}{\beta_0}$` l'ordonnée à l'origine, 
- `$\class{rouge}{\beta_l}$` l'effet de la variable `$x{,l}$` sur la variable à expliquer,
- `$\class{rouge}{\sigma^2}$` la variance.

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

- `$p+1$` paramètres de moyenne  `$(\class{rouge}{\beta_0}, \class{rouge}{\beta_1}, \ldots, \class{rouge}{\beta_p})$`; 
- 1 paramètre de variance `$\class{rouge}{\sigma^2}$`

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

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

Sauf que le modèle complet est bien plus complet ....

--
<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 }1\leq l\leq p, \class{rouge}{\beta_l} =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 }l,\  1\leq l\leq p, \class{rouge}{\beta_l} \ne 0  \right\rbrace\\
    & =\left \lbrace  M_{comp} \mbox{ est préférable à } M0 \right\rbrace.
\end{align}$$`

Sous `$H_0$`,

`$$F= \frac{ \frac{\class{alea}{RSS_0 - RSS}}{p} }{\frac{\class{alea}{RSS}}{n-(p+1)}}  \underset{H_0}{\sim}\mathcal{F}(p, n-(p+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
# Test de l'effet des variables

---
template: test_variable

## Test sur les paramètres

Tester la nullité du paramètre `$\class{rouge}{\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}$`.

`$$\begin{align} 
H_0 & =\left \lbrace \mbox{Dans le modèle complet la variable } l \mbox{ n'a pas d'influence}\right\rbrace\\
    & =\left \lbrace  \mbox{Dans le modèle complet la variable } \class{rouge}{\beta_l=0} \right\rbrace\\
    & =\left \lbrace  M_{comp} \mbox{ est équivalent à } M_{-l} \right\rbrace,
\end{align}$$`

avec 
`$$M_{-l}: \class{alea}{Y_k} = \class{rouge}{\beta_0} + \class{rouge}{\beta_1} x_{k,1} + \ldots + \class{rouge}{\beta_{l-1}} x_{k,l-1} + \class{rouge}{\beta_{l+1}} x_{k,l+1} + \ldots + \class{rouge}{\beta_{p}} x_{k,p} + \class{alea}{E_{k}}$$`

---
count: false
template: test_variable

## Equivalence des tests sur l'exemple de la pollution : effet de la population
.pull-left[

```r
Mcomp <- lm(SO2 ~ temp + manuf + pop +wind + precip + days, data = usdata)
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
```
]
--
.pull-right[
### Lien entre les statistiques  de tests

```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
```
]
---
template: test_variable
## Vigilance sur l'interprétation des tests et des estimations

--
<img src="regression_files/figure-html/ggpairs-1.png" width="60%" height="30%" style="display: block; margin: auto;" />

---
name: ajust
# Ajustement et choix de modèles

## Un critère pour mesurer l'ajustement du modèle

### Le coefficient de détermination `$R^2$` pour un modèle

`$$\class{fixe}{R^2(M)}= 1 -\frac{\class{fixe}{RSS_M}}{\class{fixe}{RSS_0}}$$`
représente la proportion de variabilité totale capturée par le modèle.
--

### Les limites de `$R^2$` pour la sélection

.pull-left[
<table>
 <thead>
  <tr>
   <th style="text-align:left;"> Model </th>
   <th style="text-align:right;"> R2 </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> M0 </td>
   <td style="text-align:right;"> 0.0000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M1 </td>
   <td style="text-align:right;"> 0.0075 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M2 </td>
   <td style="text-align:right;"> 0.8639 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M3 </td>
   <td style="text-align:right;"> 0.8702 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M4 </td>
   <td style="text-align:right;"> 0.8703 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M5 </td>
   <td style="text-align:right;"> 0.8767 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M6 </td>
   <td style="text-align:right;"> 0.8818 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M7 </td>
   <td style="text-align:right;"> 0.9011 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M8 </td>
   <td style="text-align:right;"> 0.9011 </td>
  </tr>
</tbody>
</table>
]

<IMG SRC="Ajustement.gif">
 
 ]

---
template: ajust

### Le coefficient de détermination `$R^2$`

`$$\class{fixe}{R^2(M)}= 1 -\frac{\class{fixe}{RSS_M}}{\class{fixe}{RSS_0}}$$`
représente la proportion de variabilité totale capturée par le modèle.

### Le critère Akaike (AIC)

Il utilise la vraisemblance `$^1$`  et lui ajoute une pénalité pour compenser le nombre de paramètre et éviter le surajustement.

Il est défini de manière générale pour les modèles statistiques mais dans le cas du modèle linéaire, il est donné par

$$ AIC(M) = p +  n \class{orange}{\ln{\frac{RSS}{n-(p+1))}}} + n + 3 + n\ln{2\pi} $$

$$ AIC(M) = p +  n \class{orange}{\ln{\hat{\sigma}^2_M}}  \class{clair}{ + n + 3 + n\ln{2\pi}} $$
.care[On va préférer le modèle avec l'AIC le plus petit]

.footnote[1: vraisemblance : une fonction qui mesure justement l'ajustement du modèle aux données]
---

## Comparaison AIC et `$R^2$`

--
.pull-left[
<table>
 <thead>
  <tr>
   <th style="text-align:left;"> Model </th>
   <th style="text-align:right;"> R2 </th>
   <th style="text-align:right;"> AIC </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> M0 </td>
   <td style="text-align:right;"> 0.0000 </td>
   <td style="text-align:right;"> 202.1468 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M1 </td>
   <td style="text-align:right;"> 0.0075 </td>
   <td style="text-align:right;"> 203.9960 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M2 </td>
   <td style="text-align:right;"> 0.8639 </td>
   <td style="text-align:right;"> 166.2655 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M3 </td>
   <td style="text-align:right;"> 0.8702 </td>
   <td style="text-align:right;"> 167.3043 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M4 </td>
   <td style="text-align:right;"> 0.8703 </td>
   <td style="text-align:right;"> 169.3007 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M5 </td>
   <td style="text-align:right;"> 0.8767 </td>
   <td style="text-align:right;"> 170.2800 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M6 </td>
   <td style="text-align:right;"> 0.8818 </td>
   <td style="text-align:right;"> 171.4397 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M7 </td>
   <td style="text-align:right;"> 0.9011 </td>
   <td style="text-align:right;"> 169.8802 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> M8 </td>
   <td style="text-align:right;"> 0.9011 </td>
   <td style="text-align:right;"> 171.8784 </td>
  </tr>
</tbody>
</table>
]

<IMG SRC="Ajustement.gif">
 
 ]

Utilisé en TD pour choisir un modèle avec un bon compromis ajustement et parcimonie (petit nombre de paramètres).

---

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

## Cas de la régression multiple

### Prédiction de la valeur moyenne pour des variables explicatives  `$(x_1, \ldots, x_p)$` particulières

`$$\hat{Y}_k = \class{fixe}{\hat{\beta}_0} +\class{fixe}{\hat{\beta}_1} x_{k,1}+ \ldots + \class{fixe}{\hat{\beta}_p} x_{k,p}$$`

---
# Ce qu'il faut absolument savoir pour les TDs

- savoir reconnaître une situtation de régression multiple
--

* On cherche à comprendre le lien entre une variable d'intérêt `$Y$` et des variables potentiellement explicatives toutes quantitatives.
  
--
- savoir écrire le modèle de régression multiple

--
- Tester l'existence de l'effet d'au moins une des variables explicatives candidates

--
- Tester l'existence de l'effet d'une variable explicative particulière

--
- Pourquoi il faut se méfier du `$R^2$`

--
- Ce que signifie prédire en utilisant le modèle et comment le faire