Régression linéaire

# Régression linéaire
## simple et multiple
### Marie-Pierre Etienne
### <a href="https://github.com/marieetienne" class="uri">https://github.com/marieetienne</a>
### 2020/09/11 (updated: 2023-10-02)

---

---
template: intro
## 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: intro
## 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
# 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

## Le modèle de régression simple

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

---
template: model
## Sous forme matricielle
 `$$\bf{Y = X\theta + E}$$`
### Forme régulière

`$$Y=\begin{pmatrix}
Y_{1}\\
Y_{2}\\
\vdots\\
Y_{k}\\
\vdots\\
Y_{n}\end{pmatrix},
 \quad
{\bf{X}} =\overset{\color{gray}{\begin{matrix}\beta_0  & \beta_1\end{matrix}}}{\begin{pmatrix}
1 & x_1\\
1 & x_2\\
\vdots & \vdots\\
1 & x_k\\
 \vdots & \vdots\\
1 & x_n\\
 \end{pmatrix}},\quad
{\bf{\theta}} =\begin{pmatrix}
\beta_0\\
\beta_1\\
\end{pmatrix}, \quad{\bf{E}} = \overset{}{\begin{pmatrix}
E_{1}\\
E_{2}\\
\vdots\\
E_{k}\\
\vdots\\
E_{n}\\
\end{pmatrix}}$$`
---
template: model

## Exercice

<p class=question> Détailler le modèle sous forme matricielle pour l'exemple de la pollution (écrire les 3 premières lignes de la matrice `$X$` )</p>

--
Attention, la réponse est dans la slide qui suit !

---
template: model

## Correction

```r
Mpop <- lm(SO2 ~ pop , data = usdata)
model.matrix(Mpop) %>% head(n = 3)
```

```
  (Intercept) pop
1           1 582
2           1 132
3           1 716
```
---
name: parametre
# Estimation des paramètres

---
template: parametre
## Estimation des paramètres du modèle version matricielle

Le modèle sous forme matricielle s'écrit

`$$\bf{Y = X\theta + E}.$$`
--

### Estimation de `$\theta$`

`$$\hat{\theta} = (X^\intercal X )^{-1} X^\intercal Y_{obs}.$$`

### Estimateur de `$\theta$`

`$$T = (X^\intercal X )^{-1} X^\intercal Y.$$`
--

### Loi de l'estimateur de `$\theta$`

`$$T  \sim \mathcal{N}_{I}\left(\theta, \sigma^2 (X^\intercal X )^{-1}\right).$$`

---
template: parametre
## Le paramètre de variance

La somme des carrés résiduelles s'écrit sous la forme

`$$RSS = || Y- X \hat{\theta} ||^2$$`

### Estimateur de la variance

$$S^2 =\frac{1}{DF_{res}} RSS, $$
est un <a class=care> estimateur sans bias de  `$\sigma^2$` </a> .

Dans le cas du modèle de régression simple  `$DF_{res}=n-2$` (n observations et 2 paramètres de moyennes à estimer, le nombre de composantes dans le vecteur `$\theta$`)

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

`$$\hat{\sigma}^2 =\frac{1}{n-2} RSS_{obs}.$$`
---
template: parametre
## Vérifier l'estimation sur l'exemple de la pollution

```r
X <- model.matrix(Mpop)
XXprimemoinsUn <- solve(t(X)%*%X)
XXprimemoinsUn %*% t(X) %*% matrix(usdata$SO2, ncol =1)
```

```
                   [,1]
(Intercept) 17.86831630
pop          0.02001359
```

```r
summary(Mpop)$coefficients
```

```
               Estimate  Std. Error  t value     Pr(>|t|)
(Intercept) 17.86831630 4.713843676 3.790604 0.0005093514
pop          0.02001359 0.005643813 3.546111 0.0010349968
```
--
## Loi de l'estimateur

```r
summary(Mpop)$sigma^2 * XXprimemoinsUn  # sigma2 Xt X
```

```
            (Intercept)           pop
(Intercept) 22.22032220 -1.938582e-02
pop         -0.01938582  3.185262e-05
```

```r
sqrt(diag(summary(Mpop)$sigma^2 * XXprimemoinsUn ))
```

```
(Intercept)         pop 
4.713843676 0.005643813 
```
---
name: prediction
# 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

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

## Intervalle de confiance pour le comportement moyen

---
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
## Sous forme matricielle
 `$$\bf{Y = X\theta + E}$$`
### Forme régulière

`$$Y=\begin{pmatrix}
Y_{1}\\
Y_{2}\\
\vdots\\
Y_{k}\\
\vdots\\
Y_{n}\end{pmatrix},
 \quad
{\bf{X}} =\overset{\color{gray}{\begin{matrix}\beta_0  && \beta_1&& \beta_2&&\ldots &&\beta_l\end{matrix}}}{\begin{pmatrix}
1 & x_1^{1} & x_1^{2} & \ldots &x_1^{l}\\
1 & x_2^{1} & x_2^{2} & \ldots &x_2^{l}\\
\vdots & \vdots& \vdots && \vdots\\
1 & x_k^{1} & x_k^{2} & \ldots &x_k^{l}\\
 \vdots & \vdots& \vdots && \vdots\\
1 & x_n^{1} & x_n^{2} & \ldots &x_n^{l}\\
 \end{pmatrix}},\quad
{\bf{\theta}} =\begin{pmatrix}
\beta_0\\
\beta_1\\
\beta_2\\
\vdots\\
\beta_l\\
\end{pmatrix}, \quad{\bf{E}} = \overset{}{\begin{pmatrix}
E_{1}\\
E_{2}\\
\vdots\\
E_{k}\\
\vdots\\
E_{n}\\
\end{pmatrix}}$$`

---
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
# Estimation des paramètres

---
template: parametre
## Estimation des paramètres du modèle version matricielle

Le modèle sous forme matricielle s'écrit

`$$\bf{Y = X\theta + E}.$$`
--

### Estimation de `$\theta$`

`$$\hat{\theta} = (X^\intercal X )^{-1} X^\intercal Y_{obs}.$$`

### Estimateur de `$\theta$`

`$$T = (X^\intercal X )^{-1} X^\intercal Y.$$`
--

### Loi de l'estimateur de `$\theta$`

`$$T  \sim \mathcal{N}_{I}\left(\theta, \sigma^2 (X^\intercal X )^{-1}\right).$$`

---
template: parametre
## Le paramètre de variance

La somme des carrés résiduelles s'écrit sous la forme

`$$RSS = || Y- X \hat{\theta} ||^2$$`

### Estimateur de la variance

$$S^2 =\frac{1}{DF_{res}} RSS, $$
est un <a class=care> estimateur sans bias de  `$\sigma^2$` </a> .

Dans le cas du modèle de régression simple  `$DF_{res}=n-l-1$` (n observations et 2 paramètres de moyennes à estimer, le nombre de composantes dans le vecteur `$\theta$`)

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

`$$\hat{\sigma}^2 =\frac{1}{n-l-1} RSS_{obs}.$$`
---
template: parametre
## Vérifier l'estimation sur l'exemple de la pollution

### Estimation

```r
X <- model.matrix(Mcomp)
XXprimemoinsUn <- solve(t(X)%*%X)
XXprimemoinsUn %*% t(X) %*% matrix(usdata$SO2, ncol =1)
```

```
                    [,1]
(Intercept) 111.72848064
temp         -1.26794109
manuf         0.06491817
pop          -0.03927674
wind         -3.18136579
precip        0.51235896
days         -0.05205019
```

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

### Loi de l'estimateur

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

```r
sqrt(diag(summary(Mcomp)$sigma^2 * XXprimemoinsUn ))
```

```
(Intercept)        temp       manuf         pop        wind      precip 
47.31810073  0.62117952  0.01574825  0.01513274  1.81501910  0.36275507 
       days 
 0.16201386 
```

---
name: modcomp
# 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
# 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

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

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

---
#Diagnostics

```r
library(ggfortify)
autoplot(Mcomp)
```

---
# Levier

## Mathématiquement

Le levier `$lev_i$` pour l'observation `$i$` est défini par
`$$H=X (X^\intercal X)^{-1}X^\intercal; \quad lev_{i}= H_{ii},$$`
--
### D'où viennt cette matrice `$H$` ?

$$ \hat{Y} = H Y$$
--
### Sur l'exemple de la régression simple

```r
h.us.lm <- hatvalues(Mpop)
h.us.lm
```

```
         1          2          3          4          5          6          7 
0.02444303 0.04132345 0.02524994 0.02504346 0.03952636 0.04521998 0.02603168 
         8          9         10         11         12         13         14 
0.02486268 0.02997079 0.02531882 0.59239806 0.02579734 0.03677543 0.03258749 
        15         16         17         18         19         20         21 
0.02440841 0.02896059 0.03093873 0.08536149 0.02575668 0.02515988 0.02440361 
        22         23         24         25         26         27         28 
0.02949202 0.03430016 0.04247944 0.02597074 0.02619528 0.02590162 0.02474115 
        29         30         31         32         33         34         35 
0.15851954 0.02497554 0.03814844 0.02440790 0.02631315 0.02852063 0.05345225 
        36         37         38         39         40         41 
0.03834126 0.03112651 0.03153591 0.02483924 0.04593530 0.02526602 
```

---

---
count: false

---
# Construire le modèle d'analyse de la covariance

## A partir des données chauve souris (bats)

Les différentes espèces de chauve souris ont des tailles de cerveau très variables, ce qui conduit à des volumes variables de la partie auditive.

Quel modèle pouvez vous proposer pour étudier l’influence du régime sur la part du cerveau dédiée à l’audition, compte tenu de la taille total du cerveau ?

Quel test pourrait permettre d’étudier cette influence ?

---
# Pause

<br><br><br><br>
<p style="color:#B40F20;font-size:35px;text-align:center;">Prenons une petite pause !!!

Correction en TP </p>