Régression linéaire simple, table NOTES

Présentation des variables

notes <- read.csv("~/Documents/ISIFAR/M1 ISIFAR/Semestre 2/Économétrie/Dataset/notes.csv", header = T, sep = ";")
summary(notes)
##     ID               ANGLAIS           ECO             MATH      
##  \xc9DITH  :  1   Min.   : 5.65   Min.   : 1.50   Min.   : 1.82  
##  \xc9DOUARD:  1   1st Qu.:10.39   1st Qu.: 6.50   1st Qu.: 7.92  
##  \xc9LIANE :  1   Median :11.87   Median : 8.51   Median :10.23  
##  \xc9LISE  :  1   Mean   :11.98   Mean   : 8.71   Mean   :10.01  
##  \xc9LODIE :  1   3rd Qu.:13.50   3rd Qu.:11.14   3rd Qu.:12.01  
##  \xc9LOISE :  1   Max.   :16.99   Max.   :16.75   Max.   :19.50  
##  (Other)   :119
Y = notes$MATH
X1 = notes$ECO
X2 = notes$ANGLAIS
YY = c(Y,Y)
X1X2 = c(X1,X2)
plot(Y,X1, lwd = 2, main = "Note de maths en fonction de la note d'éco", xlab = "Note économie", ylab = "Note maths", col = "blue4")

plot(Y,X2, lwd = 2, main = "Note de maths en fonction de la note d'anglais", xlab = "Note anglais", ylab = "Note maths", col = "red4")

lm 1 et régression linéaire 1:

Nous allons donc proposer un modèle de régression linéaire simple tq :

\[ Y_i = a.1 + b.X1_i + u_i \] où \(u\) est l’erreur tel que :

\[ \begin{aligned} &\mathbb E[U_i] =0 \\ &Var(U_i) = \sigma^2 \\ &\mathbb E[U_iU_j] = 0; \; \forall i \neq j \\ \end{aligned} \]

reg1 = lm(Y~X1)
summary(reg1)
## 
## Call:
## lm(formula = Y ~ X1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7647 -1.4385  0.2328  1.9853  7.4215 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.56599    0.66432   6.873 2.78e-10 ***
## X1           0.62500    0.07174   8.712 1.68e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.522 on 123 degrees of freedom
## Multiple R-squared:  0.3816, Adjusted R-squared:  0.3765 
## F-statistic: 75.89 on 1 and 123 DF,  p-value: 1.675e-14
plot(X1,Y, lwd = 2, main = "Note de maths en fonction de la note d'éco", xlab = "Note économie", ylab = "Note maths", col = "blue4")

abline(reg1, col ='black', lwd = 2)

legend("bottomright", c("Maths ~ Eco", "Y ajusté"), lty = c(3,1), lwd = c(2,2), col = c("blue4", "black"))

Pb. de test : la variable est-elle significative ?

\(H_0=\) : n’est pas signficative

\(H_1=\) : est signficative

Notre règle de décision sera de rejetter \(H_0\) si \(\alpha > p-value\). La p-value étant très proche de \(0\), on rejette \(H_0\). La variable est donc significative.

Pb. de test : la régression est-elle utile ?

\(H_0\) : la régression est inutile.

\(H_1\) : la régression utile.

ou en termes mathématiques : si \(H_0\) est vraie alors \(\mathcal F\) est toute petite et si \(H_1\) est vraie alors \(\mathcal F\) est grande. Ici la statistique de test est \(\mathcal F\), tq :

\[ \mathcal F = \frac{R^2/1}{(1-R^2)/123} \]

Ici la régression est donc utile.

Intervalle de confiance de \(\widehat b\) est donné par:

qt(0.95,123)
## [1] 1.657336
reg1$coefficients[2] + qt(0.975,123) * 0.07
##        X1 
## 0.7635612
reg1$coefficients[2] + qt(0.025,123) * 0.07
##        X1 
## 0.4864398

Verifions la gaussianité des résidus pour valider l’hypothèse sur les erreurs

shapiro.test(reg1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  reg1$residuals
## W = 0.97858, p-value = 0.04446

Pour un seuil \(\alpha = 0.05%\) la p-value est inférieure, l’hypothèse de gaussianité des résidus est alors rejetée et notre modélisation tombe à l’eau.

plot(reg1$residuals, main = "Graphe des residus",lwd = 2)

lm 2 et régression linéaire 2:

Nous allons donc proposer un modèle de régression linéaire simple tq :

\[ Y_i = a.1 + b.X2_i + u \] où \(u\) est l’erreur tel que :

\[ \begin{aligned} &\mathbb E[u_i] =0 \\ &Var(u_i) = \sigma^2 \\ &\mathbb E[u_i.u_j] = 0; \; \forall i \neq j \\ \end{aligned} \]

reg2 = lm(Y~X2)
summary(reg2)
## 
## Call:
## lm(formula = Y ~ X2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.3174 -2.0557  0.2523  2.0466  9.5474 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.3512     1.6766   6.770 4.68e-10 ***
## X2           -0.1120     0.1379  -0.812    0.418    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.198 on 123 degrees of freedom
## Multiple R-squared:  0.005334,   Adjusted R-squared:  -0.002752 
## F-statistic: 0.6597 on 1 and 123 DF,  p-value: 0.4182
plot(X2,Y, lwd = 2, main = "Note de maths en fonction de la note d'éco", xlab = "Note économie", ylab = "Note maths", col = "red4")
abline(reg2, lwd = 2, col = 'black')
legend("bottomright", c("Maths ~ Anglais", "Y"))

Prédiction

newdata = data.frame(X1 = 12)
pred1= predict.lm(reg1,newdata, level = 0.95, interval = "prediction") 
pred1
##        fit      lwr      upr
## 1 12.06599 7.032375 17.09961
newdata2 = data.frame(X2 = 6)
pred2= predict.lm(reg2,newdata2, level = 0.95, interval = "prediction") 
pred2
##        fit      lwr      upr
## 1 10.67933 4.116845 17.24182

La prévision effectuée avec la première mdoélisation est effectivement plus précise avec un intervalle de confiance de niveau \(95%\) égal à : \([7.032,17.099]\). L’intervalle est tout de même large pour l’échelle des variables sur lesquels nous travaillons, à savoir des notes d’étudiants.

Table DATA

Présentation du jeu de données :

data <- read.csv("~/Documents/ISIFAR/M1 ISIFAR/Semestre 2/Économétrie/Dataset/data.csv",header = T, sep = ";")
summary(data)
##        Y                X        
##  Min.   :0.3203   Min.   :1.101  
##  1st Qu.:1.8267   1st Qu.:3.194  
##  Median :5.7433   Median :4.660  
##  Mean   :4.7895   Mean   :4.276  
##  3rd Qu.:6.9496   3rd Qu.:5.351  
##  Max.   :9.9676   Max.   :6.184
X = data$X
Y = data$Y
plot(X,Y, lwd = 2, col = "black")

A première vue l’existence d’un lien linéaire ici semble être également une bonne piste quand on regarde le nuage de points. Modélisons alors notre variable cible \(Y\) par une régression linéaire simple, à savoir nous poserons :

\[ Y_i = a.1 + b.X_i + u_i \] où \(u\) est l’erreur tel que :

\[ \begin{aligned} &\mathbb E[u_i] =0 \\ &Var(U_i) = \sigma^2 \\ &\mathbb E[u_i.u_j] = 0; \; \forall i \neq j \\ \end{aligned} \]

MRLS

reg = lm(Y~X)
summary(reg)
## 
## Call:
## lm(formula = Y ~ X)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1149 -0.5066 -0.1983  0.8817  1.9260 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3.1737     0.7115  -4.461 0.000151 ***
## X             1.8622     0.1571  11.855 9.32e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.218 on 25 degrees of freedom
## Multiple R-squared:  0.849,  Adjusted R-squared:  0.8429 
## F-statistic: 140.5 on 1 and 25 DF,  p-value: 9.317e-12
plot(X,Y, col ='black', lwd = 2, main = "Nuage de points et régression linéaire")
abline(reg,col = 'red', lwd = 2)
legend("topleft",c("Y~X", "Droite de régression"), col = c('black', 'red'), lwd = c(3,2), lty = c(3,1))

Pb. de test : la variable est-elle significative ?

\(H_0=\) : n’est pas signficatif

\(H_1=\) : est signficatif

Notre règle de décision sera de rejetter \(H_0\) si \(\alpha > p-value\). La p-value étant très proche de \(0\), ici de l’ordre de \(10^{-12}\), on rejette \(H_0\). La variable est donc significative. De plus \(\mathcal R^2 = 0.85\), la variable explicative explique donc \(85%\) de la variable cible, ou encore, la valeur du \(cos^2 \theta\) est proche de \(1\), impliquant une forte corrélation de nos variables.

Pb de test : la régression est-elle utile ?

\(H_0\) : la régression est inutile.

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

Si l’on se pose une règle de décision sur notre variable statistique \(\mathcal F = 140\), on la remarque très éloignée de 0 ce qui nous permet d’affirmer que la régression a du sens. Esperons que l’analyse de nos résidus confirme nos hyptohèses de sur les erreurs.

Résidus : suivent-ils une loi gaussienne ?

Effectuons un shapiro.test, avec un seuil \(\alpha = 5%\) où les hypothèses sont :

\(H_0\) : sont gaussiens

\(H_1\) : ne sont pas gaussiens

shapiro.test(reg$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  reg$residuals
## W = 0.93007, p-value = 0.06936

La p-value est au-dessus de notre seuil fixé et notre indice \(W\) relativement proche de 1 ce qui est rassurant quant à la gaussianité des résidus car nous validerons alors ce test.

Régression sur un échantillon cible plus restreint :

Y8 = Y[Y<=8]
X8 = X[X<=6.02]
reg8 = lm(Y8~X8)
summary(reg8)
## 
## Call:
## lm(formula = Y8 ~ X8)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.761 -1.323 -0.218  1.069  6.458 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.3846     1.4401   0.267   0.7920  
## X8            0.9131     0.3433   2.660   0.0146 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.288 on 21 degrees of freedom
## Multiple R-squared:  0.252,  Adjusted R-squared:  0.2164 
## F-statistic: 7.076 on 1 and 21 DF,  p-value: 0.01465
plot(X8,Y8, col ='black', lwd = 2, main = "Nuage de points et regression lineaire sur echantillon restreint")
abline(reg8,col = 'blue', lwd = 2)
legend("topleft",c("Y8~X8", "Droite de regression"), col = c('black', 'blue'), lwd = c(2,2), lty = c(3,1))

Pb. de test : la variable est-elle significative ?

\(H_0=\) : n’est pas signficatif

\(H_1=\) : est signficatif

Notre règle de décision sera de rejetter \(H_0\) si \(\alpha > p-value\). La p-value étant très proche de \(0\), ici de l’ordre de \(10^{-1}\), on rejette \(H_0\). La variable est donc significative. De plus \(\mathcal R^2 = 0.25\), la variable explicative explique donc \(25%\) de la variable cible, ou encore, la valeur du \(cos^2 \theta\) est de \(0.25\), impliquant une faible corrélation de nos variables.

Pb de test : la régression est-elle utile ?

\(H_0\) : la régression est inutile.

\(H_1\) : la régression utile.

Si l’on se pose une règle de décision sur notre variable statistiques \(\mathcal F = 7\), on la remarque un peu éloignée de 0 ce qui nous permet d’affirmer que la regression a du sens.

Comparaison des régressions :

hatY = reg$coefficients[1] + reg$coefficients[2]*X
hatY8 = reg8$coefficients[1] + reg$coefficients[2]*X8
reg8$coefficients[1]
## (Intercept) 
##   0.3845905

Les droites de regression 1 et 2 sont données par :

\[ \widehat Y_{1,i}= -3.173746 + 1.862178 X_i\\ \widehat Y_{2,i} = 0.3845905 + 0.9121091 \times X_i \]

Que nous pouvons ainsi comparer en les intégrant dans le nuage de points :

plot(X,Y, col ='black', lwd = 2, main ="Tout")
abline(reg, col = "red", lwd = 2)
abline(reg8, col = 'blue', lwd =2)
legend('topleft', c('Observations',"Y chapeau 1", "Y chapeau 2"), col = c("black", "red", "blue"), lwd =c(3,2,2), lty = c(3,1,1))

Prédiction :

Calcul à la main à faire \

pred= predict.lm(reg,data.frame(X = 7), level = 0.95, interval = "prediction") 
pred
##        fit      lwr      upr
## 1 9.861498 7.158424 12.56457
pred8= predict.lm(reg8,data.frame(X8 = 7), level = 0.95, interval = "prediction") 
pred8
##        fit      lwr      upr
## 1 6.776354 1.452098 12.10061

Pour la première régression l’intervalle de confiance pour un même niveau est plus petit, le modèle 1 semble être meilleur pour effectuer des prévisions, le modèle 2 présente, pour un intervalle de confiance de niveau \(\alpha = 5%\) une plage de plus de \(10\) ce qui est relativement conséquent pour l’échelle de nos variables.

Graphe des résidus en fonction de Y

plot(reg$residuals,Y, main= "Residus en fonction de la cible", lwd =3)

plot(reg8$residuals,Y8, main= "Residus en fonction de la cible", lwd =3)

La gaussianité avait déjà été prouvée pour la première modélisation. Effectuons un shapiro.test pour l’ensemble restreint :

shapiro.test(reg8$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  reg8$residuals
## W = 0.94994, p-value = 0.2918
mean(reg8$residuals)
## [1] -1.571622e-17

On remarque que les résidus sont gaussiens.

Intervalle de confiance et prédictions :

pred1 = predict.lm(reg, level = 0.95, interval = 'prediction')
## Warning in predict.lm(reg, level = 0.95, interval = "prediction"): predictions on current data refer to _future_ responses
pred2 = predict.lm(reg, level = 0.95, interval = 'confidence')
pred3 = predict.lm(reg, level = 0.95, interval = 'none')
matplot(pred1, type = "l")

matplot(pred2, type ="l")

matplot(pred3, type ="l")

pred1;pred2;pred3
##           fit         lwr       upr
## 1   6.9022595  4.32060855  9.483910
## 2   6.2758192  3.70738067  8.844258
## 3   4.6805474  2.12505066  7.236044
## 4   2.9079124  0.33166306  5.484162
## 5   4.6260731  2.07048873  7.181657
## 6   7.8700277  5.25916447 10.480891
## 7   8.0416776  5.42453923 10.658816
## 8   5.8700395  3.30772778  8.432351
## 9   4.3891221  1.83274895  6.945495
## 10  8.2496847  5.62450811 10.874861
## 11  0.2462161 -2.42832950  2.920762
## 12  1.6640019 -0.94847624  4.276480
## 13  6.0970166  3.53151471  8.662518
## 14  8.3421735  5.71327111 10.971076
## 15  0.1777662 -2.50031291  2.855845
## 16  0.6930159 -1.95965957  3.345691
## 17  4.5408375  1.98504570  7.096629
## 18  6.6782991  4.10189190  9.254706
## 19  1.2242468 -1.40517226  3.853666
## 20  4.3056575  1.74884863  6.862466
## 21  8.1475424  5.52637217 10.768713
## 22 -1.1232815 -3.87743775  1.630875
## 23  6.4046297  3.83384514  8.975414
## 24  5.9577257  3.39425299  8.521198
## 25  8.0021454  5.38648108 10.617810
## 26  5.5045331  2.94608936  8.062977
## 27  2.6417067  0.05918176  5.224232
##           fit        lwr        upr
## 1   6.9022595  6.2956804 7.50883857
## 2   6.2758192  5.7281943 6.82344401
## 3   4.6805474  4.1972461 5.16384860
## 4   2.9079124  2.3247507 3.49107415
## 5   4.6260731  4.1423087 5.10983750
## 6   7.8700277  7.1491787 8.59087663
## 7   8.0416776  7.2984213 8.78493395
## 8   5.8700395  5.3519103 6.38816861
## 9   4.3891221  3.9012077 4.87703658
## 10  8.2496847  7.4786018 9.02076770
## 11  0.2462161 -0.6791006 1.17153289
## 12  1.6640019  0.9373254 2.39067849
## 13  6.0970166  5.5633345 6.63069863
## 14  8.3421735  7.5584998 9.12584725
## 15  0.1777662 -0.7577147 1.11324710
## 16  0.6930159 -0.1670445 1.55307626
## 17  4.5408375  4.0559784 5.02569652
## 18  6.6782991  6.0944407 7.26215750
## 19  1.2242468  0.4388417 2.00965193
## 20  4.3056575  3.8154656 4.79584933
## 21  8.1475424  7.3902115 8.90487336
## 22 -1.1232815 -2.2583557 0.01179266
## 23  6.4046297  5.8461048 6.96315459
## 24  5.9577257  5.4338851 6.48156623
## 25  8.0021454  7.2640962 8.74019471
## 26  5.5045331  5.0058840 6.00318217
## 27  2.6417067  2.0314185 3.25199487
##          1          2          3          4          5          6 
##  6.9022595  6.2758192  4.6805474  2.9079124  4.6260731  7.8700277 
##          7          8          9         10         11         12 
##  8.0416776  5.8700395  4.3891221  8.2496847  0.2462161  1.6640019 
##         13         14         15         16         17         18 
##  6.0970166  8.3421735  0.1777662  0.6930159  4.5408375  6.6782991 
##         19         20         21         22         23         24 
##  1.2242468  4.3056575  8.1475424 -1.1232815  6.4046297  5.9577257 
##         25         26         27 
##  8.0021454  5.5045331  2.6417067
matplot(X,pred1, type= "l")