Comment choisir une modélisation linéaire. À la main :

Cible = Poids ; Variable explicative = Sport

Présentation du jeu de données :

load("~/Documents/ISIFAR/M1 ISIFAR/Semestre 2/Économétrie/Dataset/physio.Rdata")
str(physio)
## 'data.frame':    20 obs. of  6 variables:
##  $ AGE   : num  2 4 6 8 10 24 48 12 34 20 ...
##  $ POIDS : num  10 15 19 25 45 48 49 50 55 58 ...
##  $ TAILLE: num  109 121 115 126 146 147 149 151 158 151 ...
##  $ SPORT : num  4.5 4.4 4.2 3.6 2.5 1.9 2.1 1.9 1.7 1.5 ...
##  $ MCDO  : num  1 2 1 2 3 2 2 3 4 6 ...
##  $ TV    : num  0.24 1.08 1.68 2.4 3.7 ...
summary(physio)
##       AGE           POIDS           TAILLE          SPORT       
##  Min.   : 2.0   Min.   :10.00   Min.   :109.0   Min.   :-0.100  
##  1st Qu.:11.5   1st Qu.:47.25   1st Qu.:146.8   1st Qu.: 1.000  
##  Median :28.0   Median :58.00   Median :151.0   Median : 1.550  
##  Mean   :27.4   Mean   :51.65   Mean   :151.6   Mean   : 1.895  
##  3rd Qu.:40.5   3rd Qu.:62.50   3rd Qu.:162.5   3rd Qu.: 2.200  
##  Max.   :58.0   Max.   :81.00   Max.   :187.0   Max.   : 4.500  
##       MCDO             TV       
##  Min.   : 1.00   Min.   :0.240  
##  1st Qu.: 2.00   1st Qu.:3.646  
##  Median : 5.00   Median :3.953  
##  Mean   : 8.40   Mean   :3.574  
##  3rd Qu.: 9.25   3rd Qu.:4.186  
##  Max.   :50.00   Max.   :5.100
Y = physio$POIDS
X = physio$SPORT

Nous avons donc \(20\) observations de \(6\) variables, on compte : AGE, POIDS, TAILLE (en entier). SPORT, MCDO, TV (en réel)

Pour la suite nous prendrons notre variable cible \(Y\) : la variable du POIDS.

Soit \(X\) = SPORT, on a alors :

plot(X,Y, main = "Poids en fonction du sport", ylab = "Poids (kg)", xlab = "Intensité sportive", lwd = 2)
reg = lm(Y ~ X)
abline(reg, col = "red", lwd = 2)
legend("topright", c("Poids~Sport","Droite de regression"), col = c("black","red"), lwd = 2)

summary(reg)
## 
## Call:
## lm(formula = Y ~ X)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1021 -1.0430  0.1494  1.1024  4.8979 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  80.1324     0.8610   93.07   <2e-16 ***
## X           -15.0303     0.3765  -39.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.155 on 18 degrees of freedom
## Multiple R-squared:  0.9888, Adjusted R-squared:  0.9882 
## F-statistic:  1594 on 1 and 18 DF,  p-value: < 2.2e-16

Les variables semblent avoir une relation linéaire avec un coefficient négatif. Nous pouvons donc essayer un modèle de régression linéaire simple. Rappelons les hypothèses :

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

\[ \begin{aligned} &\mathbb E[U] = 0 \; (centrées)\\ & Var(U_i) = \sigma^2 \; (homoscédastiques)\\ & \mathbb E[U_iU_j] = 0 \; \forall i ≠ j \; (décorrélées)\\ \end{aligned} \]

\(a\), \(b\) et \(\sigma^2\) sont les paramètres réels du modèle.

Il faut donc déterminer a et b :

\[ \left \{ \begin{array}{c @{=} c} 20a + (\sum_{i=1}^{20} X_i)b = \sum_{t=1}^{20}Y_t \\ (\sum_{i=1}^{20}X_i)a + (\sum_{i=1}^{20}X_i^2)b = \sum_{i=1}^{20} X_iY_i \end{array} \right. \Leftrightarrow \left(\begin{array}{cc} 20 & \sum_{i=1}^{20} X_i\\ \sum_{t=i}^{20}X_i & \sum_{i=1}^{20}X_i^2 \end{array}\right) \left(\begin{array}{cc} \widehat a\\ \widehat b \end{array}\right) = \left(\begin{array}{cc} \sum_{i=1}^{20}Y_i\\ \sum_{i=1}^{20} X_iY_i \end{array}\right) \Leftrightarrow \left(\begin{array}{cc} 20 & 37,9\\ 37,9 & 104,59 \end{array}\right) \left(\begin{array}{cc} \widehat a\\ \widehat b \end{array}\right) = \left(\begin{array}{cc} 1033\\ 1465 \end{array}\right) \]

M = matrix(c(20,37.9,37.9, 104.59),2,2)
b = c(sum(Y), sum(Y*X))
solve(M,b)
## [1]  80.13239 -15.03029
a_hat = 80.13239
b_hat = -15.03039

Nous avons donc avons \(\widehat a = 80.13239\) et \(\widehat b = -15.03029\)

residu = Y - a_hat - b_hat*X
mean(residu)
## [1] 0.00019905

Vérifions l’égalité de la variance

\[ \frac 1n \sum^n_{i=1}(Y_i - \bar Y)^2 = \frac 1n \sum^n_{i=1}(\widehat Y_i - \bar{\widehat Y}) + \frac 1n\sum^n_{i=1}(\widehat u - \bar{\widehat u}) \]

Y_hat = a_hat + b_hat *X
var(Y) - var(residu) - var(Y_hat)
## [1] -0.005324088

Cible : Poids ; Variable explicative : Taille

SCE<- function(X,Y){
n = length(Y)
plot(X,Y, main = "Poids en fonction de la taille", ylab = "Poids", xlab = "Taille", lwd = 2)
M = matrix(c(n,sum(X), sum(X), sum(X^2)),2,2)  
b = c(sum(Y), sum(X*Y))  
solve(M,b)
a_hat = solve(M,b)[1]
b_hat= solve(M,b)[2]
Y_hat = a_hat + b_hat*X
abline(c(a_hat,b_hat), col = 'red', lwd = 2)
legend('bottomright', c('cible~regresseur','Droite de regression'), col = c('black','red'),lty = c(1,3), lwd = 2)
Res = Y - a_hat - b_hat*X
print("Ainsi la part de variance expliquée de notre cible par la variable explicative est égale à :")
return (var(Y_hat)/var(Y))
}
SCE(physio$TAILLE,physio$POIDS)

## [1] "Ainsi la part de variance expliquée de notre cible par la variable explicative est égale à :"
## [1] 0.9414238

Ainsi la part de variance expliquée de notre cible par la variable explicative st égale à :

Cible : Poids ; Variable explicative : McDo

SCE(physio$MCDO,physio$POIDS)

## [1] "Ainsi la part de variance expliquée de notre cible par la variable explicative est égale à :"
## [1] 0.4042035

Cible : Poids ; Variable explicative : Age

SCE(physio$AGE, physio$POIDS)

## [1] "Ainsi la part de variance expliquée de notre cible par la variable explicative est égale à :"
## [1] 0.5718709

Cible : Poids ; Variable explicative : TV

SCE(physio$TV, physio$POIDS)

## [1] "Ainsi la part de variance expliquée de notre cible par la variable explicative est égale à :"
## [1] 0.9330779

Remarque interessante : POIDS/TV

Comment choisir une modélisation linéaire procédure R

Evaluation

load("~/Documents/ISIFAR/M1 ISIFAR/Semestre 2/Économétrie/Dataset/Models.Rdata")
X = Models[[1]]
Y1 = Models[[2]]
Y2 = Models[[3]]

XX = c(X,X)
YY = c(Y1,Y2)
coul = c(rep("blue",18),rep("red",18))
plot(XX,YY, col = coul, lwd = 2)
legend("topleft", c("Y1","Y2"), lwd = 2,lty = 2, col = c("blue","red"))

Le lien entre la variable \(X\) et \(Y1\) ne semble pas linéaire. Contrairement au lien entre \(X\) et \(Y2\) qui semble linéaire. Nous pourrons proposer un modèle de regression lineaire simple.

Regression linéaire

reg2 = lm(Y2~X)
reg1 = lm(Y1~X)
sum(Y1)
## [1] 8.579167
plot(XX,YY, col = coul, lwd = 2)
abline(c(reg1[1], reg1[2]), col =  "blue4", lwd = 2)
abline(c(reg2[1],reg2[2]), col ="red4", lwd = 2)
print('Moyenne des residus Y1 & Y2')
## [1] "Moyenne des residus Y1 & Y2"
mean(reg1$residuals)
## [1] -2.890453e-18
mean(reg2$residuals)
## [1] 1.263373e-19
legend("topleft", c("Y1","Y2","Reg. lin. Y1", "Reg. lin. Y2"), lty= c(3,3,1,1), lwd = 2, col = c("blue","red", "blue4", "red4"))

Bizarrement, la régression linéaire sur une variable qui ne semblait pas linéairement liée donne également une régression linéaire quasi similaire a la variable Y2. Faisons une brève étude des résidus afin d’en savoir plus, fixons notre seuil d’erreur à 5% et effectuons un test de shapiro :

plot(reg1$residuals, col = "blue", lwd = 2)
abline(c(0,0), lwd =2)
legend("topright", "Residus Y1", col = "blue", lwd = 2)

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

Remarquons qu’avec notre seuil \(\alpha = 0,06\), les résidus ne semblent pas issues d’une population normalement distribuée. Bien que l’indice statistique \(W > 0,8\) la conclusion sera que les residus ne sont pas gaussiens, et de fait, la regression lineaire ne sera pas valable.

plot(reg2$residuals, col = "red", lwd =2)
abline(c(0,0), lwd =2)
legend("topright", "Residus Y2", col = "red",lty = 3, lwd = 2)

shapiro.test(reg2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  reg2$residuals
## W = 0.94643, p-value = 0.3718

Ici les résidus semblent bien homoscédastiques. De plus le test de Shaprio-Wilk est satisfaisant ; effecivement nous avons une \(p-value = 0,37 > 0,05\), ainsi on acceptera l’hypothèse \(H_0\) de gaussianité de nos residus.