Introduction

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

library(coursesdata)
data(usdata)

La taille d'une ville est elle liée à la pollution en SO2 ?

Introduction

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.

La variable x permet elle d'expliquer la variabilité de la variable Y ?

Le modèle de régression simple

Le modèle de régression simple

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

Le modèle de régression simple

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$

Le modèle de régression simple

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

Le modèle de régression simple

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

Le modèle de régression simple

Exercice

Détailler le modèle sous forme matricielle pour l'exemple de la pollution (écrire les 3 premières lignes de la matrice $X$ )

Le modèle de régression simple

Correction

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

  (Intercept) pop
1           1 582
2           1 132
3           1 716

Estimation des paramètres

Estimation des paramètres

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).$$

Estimation des paramètres

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 estimateur sans bias de $\sigma^2$ .

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}.$$

Estimation des paramètres

Vérifier l'estimation sur l'exemple de la pollution

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

                   [,1]
(Intercept) 17.86831630
pop          0.02001359

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

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

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

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

(Intercept)         pop 
4.713843676 0.005643813

Prediction avec un modèle de régression simple

Prediction avec un modèle de régression simple

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

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

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

       fit      lwr      upr
1 24.53284 17.28453 31.78115

Prediction avec un modèle de régression simple

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

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

       fit       lwr      upr
1 24.53284 -17.90225 66.96793

Prediction avec un modèle de régression simple

Sur l'exemple de la polution

Intervalle de confiance pour le comportement moyen

Modèle de régression multiple

Modèle de régression multiple

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

Quelles sont les variables liées à la pollution en SO2 ?

Modèle de régression multiple

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$

Modèle de régression multiple

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

Modèle de régression multiple

Sur l'exemple de la pollution

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

Estimation des paramètres

Estimation des paramètres

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).$$

Estimation des paramètres

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 estimateur sans bias de $\sigma^2$ .

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}.$$

Estimation des paramètres

Vérifier l'estimation sur l'exemple de la pollution

Estimation

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

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

Estimation des paramètres

Vérifier l'estimation sur l'exemple de la pollution

Loi de l'estimateur

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

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

Test du modèle complet

Test du modèle complet

Pollution

La pollution en SO2 dans les villes américaines est elles liées à l'une au moins des variables caractérisiques des villes ?

On va à la pêche ....

Test du modèle complet

Sous forme de comparaison de modèle

Le modèle Mcomp est il plus pertinent que le modèle M0 ?

Test du modèle complet

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

Test de l'effet des variables

Test de l'effet des variables

Equivalence des tests sur l'exemple de la pollution

                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

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

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

Test de l'effet des variables

Vigilance sur l'interprétation des tests

                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

Diagnostics

library(ggfortify)
autoplot(Mcomp)

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

Prenons une petite pause !!! Correction en TP