CC BY-NC-ND 3.0
On dispose de l’observation de d’un échantillon de \(n\) couples aléatoires (\(X_i\), \(Y_i\)). Le problème est l’étude de la loi de (\(X_i\), \(Y_i\)). La question de la dépendance entre X et Y se pose.
La régression consiste à chercher une fonction \(f\) pour que \(Y_i\) soit le plus proche possible de \(f(X_i)\).
Dans le cas de la régression linéaire simple, \(f\) est \(f(x) = \beta_2x + \beta_1\). Pour estimer \(\beta_2\) et \(\beta_1\), on utilise la méthode des moindres carrés.
chapitre basé essentiellement sur le livre de Cornillon et Matzner-Lober (2011) : Régression avec R (Springer eds.)
Nous cherchons une fonction \(f\) telle que \(y_{i}=f(x_{i})\).
Pour chaque mesure i
, nous allons chercher à minimiser la différence entre \(y_i\) et \(f(x_i)\), soit :
\[ \sum_{i = 1}^{n} l(y_i - f(x_i)) \]
\(l()\) est appelée fonction de coût.
La première étape est de définir le critère de qualité et la seconde étape de définir la (ou les) fonction à utiliser.
Le choix du critère de qualité que nous allons retenir est la fonction de coût quadratique qui correspond à la différence au carré entre une observation \((x_i, y_i)\) et le point sur la droite de régression \((x_i, f(x_i))\). Nous cherchons donc à minimiser :
\[ \sum_{i = 1}^{n} (y_i - f(x_i))^2 \]
La conséquence du choix de cette fonction de coût est que les points éloignés de la droite de régression vont avoir plus d’importance que les autres.
Dans le cas de la régression linéaire simple, la fonction utilisée est la droite d’équation \(y = \beta_1 + \beta_2x\).
On suppose que nos observations \(x\) et \(y\) sont reliées sous la forme :
\[Y=\beta_1 + \beta_2X\]
On suppose aussi que la relation entre \(x\) et \(y\) est perturbée par un bruit, et notre modèle devient :
\[Y=\beta_1 + \beta_2X + \epsilon\]
Faire une régression de \(y\) sur \(x\) correspond à étudier la dépendance entre \(y\) et \(x\).
\(\beta_1\) et \(\beta_2\) sont les paramètres de notre modèle, et \(\epsilon\) le bruit ou erreur aléatoire. Pour calculer les valeurs de \(\beta_1\) et \(\beta_2\) nous allons utiliser les observations \(x_i\) et \(y_i\) :
\[y_i=\beta_1 + \beta_2x_i + \epsilon_i\]
On va chercher à minimiser la quantité \(S(\beta_1, \beta_2)\) tel que :
\[ S(\beta_1, \beta_2) = \sum_{i = 1}^{n} (y_i - f(x_i))^2 \]
# graphique des valeurs de moindres carrés
par(mfrow = c(1, 2))
plot(myMCb1, type = 'o', axes = FALSE, ylab = "S(b1, b2)",
xlab = "estimateur de b1")
axis(1, at = seq_along(estB1), labels = estB1)
plot(myMCb2, type = 'o', axes = FALSE, ylab = "S(b1, b2)",
xlab = "estimateur de b2")
axis(1, at = seq_along(estB1), labels = estB1)
par(mfrow = c(1, 1))
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 98.26464 93.74361 89.64324 85.96354 82.70451 79.86615 77.44846
## [2,] 89.94520 85.45754 81.39055 77.74423 74.51857 71.71359 69.32927
## [3,] 82.12576 77.67147 73.63786 70.02491 66.83263 64.06102 61.71008
## [4,] 74.80631 70.38540 66.38516 62.80559 59.64669 56.90846 54.59089
## [5,] 67.98687 63.59933 59.63247 56.08628 52.96075 50.25589 47.97171
## [6,] 61.66742 57.31327 53.37978 49.86696 46.77481 44.10333 41.85252
## [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 75.45144 73.87508 72.71940 71.98438 71.67004 71.77636 72.30335
## [2,] 67.36563 65.82265 64.70034 63.99870 63.71773 63.85743 64.41779
## [3,] 59.77981 58.27021 57.18128 56.51301 56.26542 56.43849 57.03223
## [4,] 52.69400 51.21777 50.16222 49.52733 49.31311 49.51956 50.14668
## [5,] 46.10819 44.66534 43.64315 43.04164 42.86080 43.10062 43.76112
## [6,] 40.02237 38.61290 37.62409 37.05596 36.90849 37.18169 37.87556
## [,15] [,16] [,17] [,18] [,19] [,20] [,21]
## [1,] 73.25101 74.61934 76.40834 78.61800 81.24834 84.29934 87.77101
## [2,] 65.39883 66.80053 68.62291 70.86595 73.52966 76.61404 80.11909
## [3,] 58.04665 59.48173 61.33748 63.61389 66.31098 69.42874 72.96716
## [4,] 51.19446 52.66292 54.55205 56.86184 59.59230 62.74343 66.31524
## [5,] 44.84228 46.34411 48.26662 50.60979 53.37362 56.55813 60.16331
## [6,] 38.99010 40.52531 42.48119 44.85773 47.65495 50.87283 54.51138
# gradient de couleur du rouge au bleu
# on ne rentre pas dans les détails ici...
myCol <- colorRampPalette(c("blue", "white", "red"))(100)
myColRankc <- (myMCb12[-1, -1] +
myMCb12[-1, -ncol(myMCb12)] +
myMCb12[-nrow(myMCb12), -1] +
myMCb12[-nrow(myMCb12), -ncol(myMCb12)])/4
myColRankc <- cut(myColRankc, 100)
On peut observer que la fonction \(S(\beta_1, \beta_2)\) est convexe, il n’y a donc qu’une seule solution. Cette solution peut être calculée analytiquement.
\[\hat{\beta_1} = \overline{y} - \hat{\beta_2}\overline{x}\]
\[\hat{\beta_2} = \frac{\sum(x_i - \overline{x})y_i} {\sum(x_i - \overline{x})^2} \]
Note : “Régression avec R” page 10 pour plus d’information sur les calculs.
estB2 <- sum((myX - mean(myX))*myY) / sum((myX - mean(myX))^2)
estB1 <- mean(myY) - estB2 * mean(myX)
fitY <- estB1 + estB2 * myX
print(paste0("estimateur de béta 1 = ", estB1))
## [1] "estimateur de béta 1 = 0.823768464579744"
## [1] "estimateur de béta 2 = 0.497004906743087"
myCol <- colorRampPalette(c("green", "blue", "red"))(101)
colRank <- (myY - fitY)^2
colRank <- round((colRank - min(colRank)) /
(max(colRank) - min(colRank)) * 100) + 1
plot(x = myX, y = myY, col = myCol[colRank], pch = 16,
axes = FALSE, panel.first = {
grid()
axis(1)
axis(2)
})
points(x = myX, y = fitY, type = 'l')
myCol <- colorRampPalette(c("green", "blue", "red"))(101)
colRank <- (myY - fitY)^2
colRank <- round((colRank - min(colRank)) /
(max(colRank) - min(colRank)) * 100) + 1
plot(x = myX, y = myY, col = myCol[colRank], pch = 16,
axes = FALSE, panel.first = {
grid()
axis(1)
axis(2)
segments(
x0 = myX,
y0 = fitY,
x1 = myX,
y1 = myY,
col = myCol[colRank], lwd = 1
)
})
points(x = myX, y = fitY, type = 'l')
Les résidus sont calculés avec \(\hat{\epsilon_i}=y_i-\hat{y_i}\). La somme des résidus est nulle (ici très proche de zéro à cause du nombre de chiffres significatifs dans les calculs).
## [1] 4.302114e-15
La variance résiduelle est estimée par la variation des résidus :
\[s^2=\frac{1}{n-2}\sum_{i=1}^{n}{e_i}^2\]
## [1] 0.0415712
Pour une nouvelle valeur de X \(x_{n+1}\), nous pouvons calculer \(\hat{y_{n+1}} = \hat{\beta_1}+\hat{\beta_2}x_{n+1}\).
## [1] 1.320773 1.817778 2.314783 2.811788
1- la distribution de \(\epsilon\) est indépendante de \(X\).
2- \(\epsilon\) est distribué normalement, centré sur zéro et de variance constante (homoscédasticité)
1- la distribution de \(\epsilon\) est indépendante de \(X\).
2- \(\epsilon\) est distribué normalement, centré sur zéro et de variance constante
##
## Shapiro-Wilk normality test
##
## data: residus
## W = 0.99324, p-value = 0.902
2- \(\epsilon\) est distribué normalement, centré sur zéro et de variance constante
2- \(\epsilon\) est distribué normalement, centré sur zéro et de variance constante
2- \(\epsilon\) est distribué normalement, centré sur zéro et de variance constante
\[SST=\sum(y_i - \overline{y})^2\]
## [1] 24.82927
\[SSM=\sum(\hat{y_i} - \overline{y})^2\]
## [1] 20.67215
\[SSR=\sum(\hat{y_i} - {y_i})^2\]
## [1] 4.15712
## [1] 4.15712
## [1] 0.8325718
lm
lm
Dans la pratique nous utilisons la fonction lm
:
##
## Call:
## lm(formula = myY ~ myX)
##
## Coefficients:
## (Intercept) myX
## 0.8238 0.4970
lm
##
## Call:
## lm(formula = myY ~ myX)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.59376 -0.14368 -0.00841 0.13534 0.60415
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.82377 0.02065 39.89 <2e-16 ***
## myX 0.49700 0.02251 22.07 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.206 on 98 degrees of freedom
## Multiple R-squared: 0.8326, Adjusted R-squared: 0.8309
## F-statistic: 487.3 on 1 and 98 DF, p-value: < 2.2e-16
lm
Coefficients:
lm
myCol <- colorRampPalette(c("green", "blue", "red"))(101)
colRank <- (myY - fitY)^2
colRank <- round((colRank - min(colRank)) /
(max(colRank) - min(colRank)) * 100) + 1
myPrev <- predict(myModel, interval = "confidence", level = 0.95)
plot(x = myX, y = myY, col = myCol[colRank], pch = 16,
axes = FALSE, panel.first = {
grid()
axis(1)
axis(2)
})
points(x = myX, y = myPrev[,1], type = 'l', lwd = 2)
points(x = myX, y = myPrev[,2], type = 'l', lty = 2, col = 2)
points(x = myX, y = myPrev[,3], type = 'l', lty = 2, col = 2)
lm
lm
## 2.5 % 97.5 %
## (Intercept) 0.7827877 0.8647493
## myX 0.4523268 0.5416830
lm
lm
##
## Call:
## lm(formula = myYsample ~ myXsample)
##
## Coefficients:
## (Intercept) myXsample
## 0.8393 0.5158
## [1] 0.7988975 0.8014190
## attr(,"conf.level")
## [1] 0.95
## [1] 0.4989557 0.5014942
## attr(,"conf.level")
## [1] 0.95
## [1] 0.02031668
## [1] 0.0204538
##
## Call:
## lm(formula = myY ~ myX)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.88834 -0.13539 0.00053 0.13433 1.02864
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7988522 0.0006327 1262.6 <2e-16 ***
## myX 0.4998063 0.0006355 786.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2001 on 99998 degrees of freedom
## Multiple R-squared: 0.8608, Adjusted R-squared: 0.8608
## F-statistic: 6.186e+05 on 1 and 99998 DF, p-value: < 2.2e-16
## 2.5 % 97.5 %
## (Intercept) 0.7976121 0.8000922
## myX 0.4985608 0.5010519