Regression lineaire simple, table NOTES

Presentation des variables

notes =  read.csv("~/data/notes.csv", header = TRUE, sep =";" )
summary(notes)
##         ID         ANGLAIS           ECO             MATH      
##  ABSOLON :  1   Min.   : 5.65   Min.   : 1.50   Min.   : 1.82  
##  ADÉLAÏDE:  1   1st Qu.:10.39   1st Qu.: 6.50   1st Qu.: 7.92  
##  ADÈLE   :  1   Median :11.87   Median : 8.51   Median :10.23  
##  ADOLPHE :  1   Mean   :11.98   Mean   : 8.71   Mean   :10.01  
##  ADRIEN  :  1   3rd Qu.:13.50   3rd Qu.:11.14   3rd Qu.:12.01  
##  ADRIENNE:  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 = 3, main = "Note de math en fonction de la note d'éco", xlab = "Note economie", ylab = "Note maths", col = "blue4")

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

lm 1 et regression lineaire 1:

Nous allons donc proposer un modele de regression lineaire simple tq :

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

\[ \begin{aligned} &1-\mathbb E[U_i] =0 \\ & 2-Var(U_i) = \sigma^2 \\ & 3- \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 = 3, main = "Note de math en fonction de la note d'éco", xlab = "Note economie", ylab = "Note maths", col = "blue4")

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

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

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\), on rejette \(H_0\). La variable est donc significative.

Pb de test : la regression 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)/48} \]

Ici la regression est significative.

Intervalle de confiance de \(\hat 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 gaussianite des residus pour valider l’hypothese 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’hypothese de gaussianite des residus est alors rejete et notre modelisation tombe alors à l’eau…

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

lm 2 et regression lineaire 2:

Nous allons donc proposer un modele de regression lineaire simple tq :

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

\[ \begin{aligned} &1-\mathbb E[U_i] =0 \\ & 2-Var(U_i) = \sigma^2 \\ & 3- \mathbb E[U_iU_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 = 3, main = "Note de math en fonction de la note d'éco", xlab = "Note economie", ylab = "Note maths", col = "red4")
abline(reg2, lwd =2, col = 'black')
legend("bottomright", c("Maths ~ Anglais", "Y"))

Prediction

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

Presentation du jeu de données :

data = read.csv("~/data/dataexo.txt", 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 = 3, col = "black")

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

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

\[ \begin{aligned} &1-\mathbb E[U_i] =0 \\ & 2-Var(U_i) = \sigma^2 \\ & 3- \mathbb E[U_iU_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 = 3, main = "Nuage de points et regression lineaire")
abline(reg,col = 'red', lwd =2)
legend("topleft",c("Y~X", "Droite de regression"), 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 regression 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 = 140\), on l’a remarque très éloignée de 0 ce qui nous permet d’affirmer que la regression à du sens. Esperons que l’analyse de nos résidus affirment nos hyptohèses de sur les erreurs…

Residus : 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.

Regression sur un echantillon 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 = 3, 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(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^{-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 regression 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 l’a remarque un peu éloignée de 0 ce qui nous permet d’affirmer que la regression à 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 :

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

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

plot(X,Y, col ='black', lwd =3, 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 regression 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 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é 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 la résidus sont gaussiens et cela avec plus de certitudes que les premiers résidus.. personnellement j’aurais dit le contraire juste à la gueule des graphes. Bon après réflexion le graphe est tracé en fonction de \(Y \le 8\),les résidus ne suivent donc certainement pas une loi \(\mathcal N\) iid, mais admettent certainement une dépendance à la variable \(Y\)

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

matplot(pred1, type = "l")

matplot(X,pred1, type= "l")