Régression linéaire multiple

Modélisation de l’ozone à l’aide de toutes les variables explicatives

library(readr)
ozo = read.csv("~/Documents/ISIFAR/M1 ISIFAR/Semestre 2/Économétrie/Dataset/ozone.csv", header = T, sep = ";")

Objectif : établir un modèle de régression multiple expliquant la pollution.

Soit un seuil d’erreur de niveau \(\alpha = 0.05\).

Ici, notre variable cible est O3 et nos features sont les T12, T15, Ne12, N12, S12, E12, W12, Vx.

Modèle de régression multiple :

\[ Y_i = b_1.X_i^1 + b_2.X_i^2 + ... + b_{p}.X_i^{p} + u_i \\ \]

\[ \underbrace{Y}_{(n \times 1)} = \underbrace{X}_{(n\times p)}.\underbrace{b}_{(p \times 1)} + \underbrace{u}_{(n \times 1)} \\ \]

\[ \begin{pmatrix} Y_1 \\ \vdots \\ Y_n \end{pmatrix} = \begin{pmatrix} X^1_1 & \cdots & X_1^{p} \\ \vdots & \ddots & \vdots \\ X_1^{n} & \cdots & X_n^{p} \end{pmatrix} \times \begin{pmatrix} b_1 \\ \vdots \\ b_{p} \end{pmatrix} + \begin{pmatrix} u_1 \\ \vdots \\ u_{n} \end{pmatrix} \]

où \(u\) est l’erreur telle que :

\[ \left\{ \begin{aligned} &\mathbb E[u] = 0_n \; (centrées)\\\\ & Var(u) = \begin{pmatrix} Var(u_1) & & Cov(u_i,u_j) \\ & \ddots \\ & & Var(u_n) \end{pmatrix} = \underbrace{\sigma^2}_{homoscédasticité} \times \underbrace{Id_n}_{décorrélation} \end{aligned} \right. \]

Ainsi nous avons les paramètres du modèle :

\(b \in \mathbb R^{p-1}\)

\(\sigma^2>0\)

\(u \sim \mathcal N(0_n, \sigma^2.Id_n)\)

feat = data.frame(ozo$T12, ozo$T15, ozo$Ne12, ozo$N12, ozo$S12, ozo$E12, ozo$W12, ozo$Vx, ozo$maxO3v)
reg = lm(ozo$O3 ~., feat)
summary(reg)

## 
## Call:
## lm(formula = ozo$O3 ~ ., data = feat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.0532  -7.7178   0.8154   7.9974  27.9616 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 54.72777   17.27890   3.167 0.002943 ** 
## ozo.T12     -0.35184    1.57308  -0.224 0.824159    
## ozo.T15      1.49719    1.53774   0.974 0.336092    
## ozo.Ne12    -4.19219    1.06376  -3.941 0.000318 ***
## ozo.N12      1.27549    1.36319   0.936 0.355061    
## ozo.S12      3.17112    1.91078   1.660 0.104817    
## ozo.E12      0.52773    1.94274   0.272 0.787294    
## ozo.W12      2.47488    2.07201   1.194 0.239342    
## ozo.Vx       0.60770    0.48575   1.251 0.218187    
## ozo.maxO3v   0.24536    0.09649   2.543 0.014964 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.21 on 40 degrees of freedom
## Multiple R-squared:  0.7505, Adjusted R-squared:  0.6944 
## F-statistic: 13.37 on 9 and 40 DF,  p-value: 1.507e-09

var(ozo$O3)

## [1] 571.2494

Nous testerons alors chaque variable une par une :

Tests nécessaires à la validation du modèle :

Test de Fisher : La régression est-elle signficative ?

\[ \left\{ \begin{aligned} &H_0 : la \; régression \; est \; non \; significative \Leftrightarrow b_i = 0 \; \forall i\in [1,p]\\ & H_1 : la \; régression \; est \; significative \Leftrightarrow \exists i\in [1,p] \; tq \; b_i ≠ 0 \end{aligned} \right. \] Indicateurs : on regarde \(R^2\) / \(F-statistic\) / \(p-value\)

On établit la règle suivante :

Si p−value \(<α\) \(\Rightarrow\) On rejette \(H_0\) / On accepte \(H_1\) \(\Rightarrow\) La régression est significative.

Si p−value \(>α\) \(\Rightarrow\) On accepte \(H_0\) / On rejette \(H_1\) \(\Rightarrow\) La régression n’est pas significative.

Test de Student :

\[ \left\{ \begin{aligned} &H_0 : la \; variable \; X^j \; est \; non \; significative \Leftrightarrow b_j = 0 \; \forall j\in [1,p]\\ & H_1 : la \; variable\; X^j \; est \; significative \Leftrightarrow \exists j\in [1,p] \; tq \; b_j ≠ 0 \end{aligned} \right. \]

Indicateurs : on regarde \(T\) / \(p\)-value

Règle de décision :

Si p−value \(<α\) \(\Rightarrow\) On rejette \(H_0\) / On accepte \(H_1\) \(\Rightarrow\) La variable est significative.

Si p−value \(>α\) \(\Rightarrow\) On accepte \(H_0\) / On rejette \(H_1\) \(\Rightarrow\) La variable n’est pas significative.

Consolider une sortie R

a/ On remarque que le sortie R nous informe qu’il y a 19 degrés de liberté, 5 variables sont impliquées dans la modélisation : nous avions alors 24 observations sur ces données.

b/ 1- Intercept

2- \(-1518.18 / 739.93 = -2.05\) Calcul de la T valeur

3- \(34.67 / 0.178 = 194,7753\) Calcul de la T valeur

4- Rien du tout

5- \(14.31 \times 2.697 = 38.59407\) Calcul de la T valeur

6- 4 (Nombre de variables explicatives)

7- 19 (Degrés de liberté)

c/ SCT = SCE + SCR \(R^2 = \frac{SCE}{SCT}\)

Ici \(R^2 = 0.86\)

Taux de remplissage des trains

La variable d’intérêt est le taux de remplissage d’un train

Audit

train = read.csv("~/Documents/ISIFAR/M1 ISIFAR/Semestre 2/Économétrie/Dataset/RemplissageTrain.csv", header = TRUE, sep =";", dec = ",")
summary(train)

##        TX             HEURE            JOUR           SEMAINE  
##  Min.   :0.1464   Min.   : 7.00   Min.   :  1.00   Min.   : 1  
##  1st Qu.:0.4680   1st Qu.:10.00   1st Qu.: 38.00   1st Qu.: 6  
##  Median :0.6246   Median :11.00   Median : 74.00   Median :11  
##  Mean   :0.6058   Mean   :12.35   Mean   : 74.03   Mean   :11  
##  3rd Qu.:0.7388   3rd Qu.:14.00   3rd Qu.:110.00   3rd Qu.:16  
##  Max.   :1.0000   Max.   :21.00   Max.   :147.00   Max.   :21  
##   TEMPERATURE      PUISSANCE      DUREESOLEIL      PLUIE       
##  Min.   :18.80   Min.   :142.0   Min.   : 18   Min.   : 0.000  
##  1st Qu.:25.60   1st Qu.:753.0   1st Qu.:560   1st Qu.: 0.000  
##  Median :27.80   Median :845.0   Median :702   Median : 0.000  
##  Mean   :27.58   Mean   :807.6   Mean   :635   Mean   : 1.183  
##  3rd Qu.:29.80   3rd Qu.:922.0   3rd Qu.:804   3rd Qu.: 0.000  
##  Max.   :35.20   Max.   :997.0   Max.   :872   Max.   :29.800

str(train)

## 'data.frame':    777 obs. of  8 variables:
##  $ TX         : num  0.146 0.325 0.342 0.273 0.24 ...
##  $ HEURE      : int  10 11 12 14 7 8 10 8 9 10 ...
##  $ JOUR       : int  1 1 1 1 2 2 2 3 3 3 ...
##  $ SEMAINE    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ TEMPERATURE: num  23.1 23.1 23.1 23.1 18.8 18.8 18.8 20 20 20 ...
##  $ PUISSANCE  : int  839 839 839 839 492 492 492 925 925 925 ...
##  $ DUREESOLEIL: int  642 642 642 642 319 319 319 745 745 745 ...
##  $ PLUIE      : num  0 0 0 0 0 0 0 0 0 0 ...

Modélisation

trainregresseur = data.frame(train$HEURE, train$JOUR, train$SEMAINE , train$TEMPERATURE, train$PUISSANCE, train$DUREESOLEIL, train$PLUIE)

reg2 = lm(train$TX ~., trainregresseur)
summary(reg2)

## 
## Call:
## lm(formula = train$TX ~ ., data = trainregresseur)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.233287 -0.041098  0.001505  0.039796  0.191818 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        1.794e-02  2.438e-02   0.736   0.4621    
## train.HEURE        6.798e-03  6.452e-04  10.536  < 2e-16 ***
## train.JOUR         1.997e-03  1.365e-03   1.463   0.1440    
## train.SEMAINE      1.170e-02  9.591e-03   1.220   0.2229    
## train.TEMPERATURE  5.946e-03  8.010e-04   7.423 3.05e-13 ***
## train.PUISSANCE    8.124e-06  2.386e-05   0.341   0.7335    
## train.DUREESOLEIL  9.134e-05  1.604e-05   5.696 1.75e-08 ***
## train.PLUIE       -9.883e-04  5.201e-04  -1.900   0.0578 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06543 on 769 degrees of freedom
## Multiple R-squared:  0.8659, Adjusted R-squared:  0.8647 
## F-statistic: 709.2 on 7 and 769 DF,  p-value: < 2.2e-16

Problème de test (Student) : les variables sont-elles significatives ?

\(H_0\) : la variable \(X_p\) n’est pas significative, ou \(b_p = 0\).

\(H_1\) : la variable \(X_p\) est significative, ou \(b_p \not = 0\).

Règle de décision : ici, on regarde l’indice statistique : \(\mathcal T = \frac{\widehat b}{\sqrt{\widehat{Var(\widehat b)}}}\ \)~\(\ \mathcal t_{1587}\). On établit la règle suivante :

t-value > \(qt_{1587}\)(\(\frac{1+ \alpha}{2}\))

qt(1.05/2,769)

## [1] 0.06272725

Problème de test (Fisher) : la régression est-elle significative ?

\(H_0\) : la régression n’est pas significative.

\(H_1\) : la régression est significative.

Règle de décision : nous regardons l’indice statistique : \(\mathcal F = \frac {\mathbb R^2/1}{(1-\mathbb R^2)/1587}\)

On a \(F = 709.2\) >

qf(0.05,1,769)

## [1] 0.003934708

Problème de test (Shapiro) : les résidus sont-ils gaussiens ?

\(H_0\) = les résidus sont gaussiens

\(H_1\) : les résidus ne sont pas gaussiens

shapiro.test(reg2$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  reg2$residuals
## W = 0.99737, p-value = 0.2491

Ainsi, la p-value est supérieure au seuil \(\alpha = 0.05\) et la statistique \(W\) est très proche de 1, on accepte alors l’hypothèse nulle et les résidus sont donc problamement issus d’une population normalement distribuée.

Prédiction

index = sample(1:777, 539, replace = FALSE)
TRAIN = train[index,]
TEST = train[-index,]

Modélisons sur notre variable TRAIN:

TRAINREG = data.frame(TRAIN$HEURE, TRAIN$JOUR, TRAIN$SEMAINE , TRAIN$TEMPERATURE, TRAIN$PUISSANCE, TRAIN$DUREESOLEIL, TRAIN$PLUIE)
regTrain = lm(TRAIN$TX~.,TRAINREG)
summary(regTrain)

## 
## Call:
## lm(formula = TRAIN$TX ~ ., data = TRAINREG)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.233352 -0.041776  0.000881  0.038652  0.190802 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        1.867e-02  2.985e-02   0.625 0.531946    
## TRAIN.HEURE        6.699e-03  7.713e-04   8.686  < 2e-16 ***
## TRAIN.JOUR         3.526e-03  1.679e-03   2.100 0.036178 *  
## TRAIN.SEMAINE      6.437e-04  1.180e-02   0.055 0.956535    
## TRAIN.TEMPERATURE  6.580e-03  9.802e-04   6.713 4.91e-11 ***
## TRAIN.PUISSANCE    8.841e-06  2.892e-05   0.306 0.759976    
## TRAIN.DUREESOLEIL  7.547e-05  1.983e-05   3.805 0.000159 ***
## TRAIN.PLUIE       -1.048e-03  5.987e-04  -1.751 0.080512 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06624 on 531 degrees of freedom
## Multiple R-squared:  0.8617, Adjusted R-squared:  0.8598 
## F-statistic: 472.5 on 7 and 531 DF,  p-value: < 2.2e-16