Suite au premier exercice sur la régression linéaire simple avec R, voici un nouvel exercice sur la régression linéaire multiple avec R.

À nouveau, je vais dans un premier temps présenter toutes les étapes comme on pourrait les faire à la main, puis je terminerai par les deux lignes de code qui permettent d’obtenir les mêmes résultats.

Le but de cet exercice est d’appliquer les formules qui permettent d’obtenir les estimateurs de paramètres de la régression, et d’effectuer les tests d’hypothèses. Il faut garder à l’esprit que lorsque l’on souhaite effectuer une régression, il ne faut pas se lancer directement dans les calculs, mais prendre son temps pour observer les données et regarder quels types de relations les lient entre-elles (ce que nous ne ferons pas dans cet exercice).

On dispose d’une variable endogène (\(y\)) dont on souhaite étudier les variations, en s’appuyant sur quatre variables exogènes (\(x_1,x_2,x_3,x_4\)).

On définit la matrice \(\boldsymbol X\) comme suit :

\( \boldsymbol X = \begin{bmatrix}
x_{11} & x_{12} & x_{13} & x_{14} & 1 \\
x_{21} & x_{22} & x_{23} & x_{24} & 1 \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
x_{n1} & x_{n2} & x_{n3} & x_{n4} & 1
\end{bmatrix}\)

Le nombre d’observations s’élève à 31.

> y <- c(5.7,5.8,6.1,6.5,6.8,6.8,7.1,21.3,18.7,14.5,7.4,9,
+        11.7,9.5,9.5,8.8,9.3,8.6,7.7,10.8,6.6,11.7,11.9,
+        10.8,7.6,11.3,10.8,9.2,11.6,12.8,12.7)
> 
> x1 <- c(11600,12490,10450,17140,14825,13730,19490,285000,
+         183900,92500,25000,22350,36600,22500,31580,28750,
+         22600,20300,19900,39800,19740,38990,50800,36200,
+         31990,47700,36950,26950,36400,50900,49300)
> 
> x2 <- c(846,993,899,1390,1195,658,1331,5474,5987,2789,
+         1597,1761,2165,1983,1984,1998,1580,1390,1396,
+         2435,1242,2972,2958,2497,1998,2496,1998,1997,
+         1984,2438,2473)
> 
> x3 <- c(32,39,29,44,33,32,55,325,300,209,74,74,101,85,
+         85,89,65,54,66,106,55,107,150,122,66,125,89,92,
+         85,97,125)
> 
> x4 <- c(650,790,730,955,895,740,1010,1690,2250,1485,1080,
+         1100,1500,1075,1155,1140,1080,1110,1140,1370,940,
+         1400,1550,1330,1300,1670,1560,1240,1635,1800,1570)
> n <- length(y)
> cste <- rep(1, n)
> 
> X <- cbind(x1 = x1, x2 = x2, x3 = x3, x4 = x4, cste =cste)
> 
> y ; X
 [1]  5.7  5.8  6.1  6.5  6.8  6.8  7.1 21.3 18.7 14.5  7.4  9.0 11.7  9.5  9.5  8.8  9.3  8.6  7.7 10.8
[21]  6.6 11.7 11.9 10.8  7.6 11.3 10.8  9.2 11.6 12.8 12.7
          x1   x2  x3   x4 cste
 [1,]  11600  846  32  650    1
 [2,]  12490  993  39  790    1
 [3,]  10450  899  29  730    1
 [4,]  17140 1390  44  955    1
 [5,]  14825 1195  33  895    1
 [6,]  13730  658  32  740    1
 [7,]  19490 1331  55 1010    1
 [8,] 285000 5474 325 1690    1
 [9,] 183900 5987 300 2250    1
[10,]  92500 2789 209 1485    1
[11,]  25000 1597  74 1080    1
[12,]  22350 1761  74 1100    1
[13,]  36600 2165 101 1500    1
[14,]  22500 1983  85 1075    1
[15,]  31580 1984  85 1155    1
[16,]  28750 1998  89 1140    1
[17,]  22600 1580  65 1080    1
[18,]  20300 1390  54 1110    1
[19,]  19900 1396  66 1140    1
[20,]  39800 2435 106 1370    1
[21,]  19740 1242  55  940    1
[22,]  38990 2972 107 1400    1
[23,]  50800 2958 150 1550    1
[24,]  36200 2497 122 1330    1
[25,]  31990 1998  66 1300    1
[26,]  47700 2496 125 1670    1
[27,]  36950 1998  89 1560    1
[28,]  26950 1997  92 1240    1
[29,]  36400 1984  85 1635    1
[30,]  50900 2438  97 1800    1
[31,]  49300 2473 125 1570    1

Le modèle que l’on estime s’écrit :
\[y_i = \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + \beta_4 x_{4i} + \beta_0 + \varepsilon_i, \quad i=1,2,\ldots, n\]
ou de manière équivalente, sous forme matricielle :
\[\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon},\]
avec \(\boldsymbol{y} = \begin{bmatrix}
y_{1} & y_{2} & \cdots & y_{n}
\end{bmatrix}^t \), \(\boldsymbol{\beta} = \begin{bmatrix} \beta_1 & \beta_2 & \beta_3 & \beta_4 & \beta_0 \end{bmatrix}^t\), \(\boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_1 & \varepsilon_2 & \ldots & \varepsilon_n \end{bmatrix}^t\) et la matrice \(\boldsymbol{X}\) définie plus haut.

Les estimateurs MCO des coefficients de la régression sont donnés par :
\[\hat{\boldsymbol\beta} = (\boldsymbol X^t \boldsymbol X)^{-1} \boldsymbol X^t \boldsymbol y.\]

Effectuons ces calculs, pas-à-pas.

> # La transposee de X
> # t(X)
> 
> # X^t X
> t(X) %*% X
               x1         x2        x3         x4    cste
x1   150316730025 4548968325 238826625 2082249275 1356425
x2     4548968325  175725810   8609270   91977255   64904
x3      238826625    8609270    436506    4366165    3010
x4     2082249275   91977255   4366165   52817950   38940
cste      1356425      64904      3010      38940      31
> 
> # Inversion de la matrice X^t X, soit (X^t X)^(-1)
> solve(t(X) %*% X)
                x1            x2            x3            x4          cste
x1    1.141294e-10 -1.880855e-09 -7.710919e-08  5.530438e-09 -5.157935e-07
x2   -1.880855e-09  4.947312e-07 -4.808522e-06 -4.210556e-07  4.228260e-05
x3   -7.710919e-08 -4.808522e-06  1.494839e-04 -2.063264e-06  1.518776e-03
x4    5.530438e-09 -4.210556e-07 -2.063264e-06  1.156298e-06 -6.125560e-04
cste -5.157935e-07  4.228260e-05  1.518776e-03 -6.125560e-04  5.882820e-01
> 
> # X^t y
> t(X) %*% y
           [,1]
x1   18756861.0
x2     762211.7
x3      37036.7
x4     421009.0
cste      308.6
> 
> # \betaChap = (X^t X)^(-1) X^t y
> (hatBeta <- (solve(t(X) %*% X)) %*% (t(X) %*% y))
              [,1]
x1    2.042054e-05
x2   -5.005933e-04
x3    2.499448e-02
x4    4.160583e-03
cste  2.456294e+00

À partir de ces coefficients, on peut calculer à présent les estimations \(\hat{\boldsymbol{y}}\), et ensuite obtenir les résidus :

> yChap <- X %*% hatBeta
> residus <- y - yChap

On peut calculer le coefficient de détermination (\(R^2\)) à l'aide de la relation suivante :
\[R^2 = \frac{SCE}{SCT},\]
avec \(SCE = \sum_{i=1}^{n}(\hat{y}_i - \bar{y})^2\) et \(SCT = \sum_{i=1}^{n}(y-\bar{y})^2\),
où \(\bar{y} = n^{-1} \sum_{i=1}^{n} y_i\) et \(\bar{y} = n^{-1} \sum_{i=1}^{n} x_i\).

> (SCE <- sum((yChap - mean(y))^2))
[1] 364.7719
> (SCT <- sum((y - mean(y))^2))
[1] 382.1368
> (R2 <- SCE / SCT)
[1] 0.9545586

La lecture du \(R^2\) nous indique que \(95.45\%\) des variations de \(y\) sont expliquées par le modèle. Il faut toutefois rester prudent. Dans cet exercice, on se précipite sur les calculs de régression, sans avoir jeté un oeil aux données, sans avoir regardé les corrélations existantes entre les variables, etc. Aussi, toutes les interprétations que je donne ici sont à prendre avec des pincettes, et donnent juste une clé de lecture dans le cas où tout va bien.

On a calculé le coefficient de détermination, calculons à présent le coefficient de corrélation ajusté, qui vient apporter une pénalité au \(R^2\), afin de prendre en compte le nombre de variables explicatives incluses dans le modèle. Il est défini comme suit :
\[R^2_a = 1 - \frac{n-1}{n-m-1}(1-R^2),\]
avec \(m\) le nombre de variables explicatives.

> m <- 4 # Nombre de variables explicatives
> (R2a <- 1 - ((n-1) / (n-m-1)) * (1-R2))
[1] 0.9475676

Afin de pouvoir effectuer des tests de significativité pour chacun des coefficients, nous avons besoin de calculer au préalable l'estimation de la variance des erreurs ainsi que les estimations de la variance des estimateurs des paramètres (les éléments diagonaux de la matrice de variance-covariance).

L'estimation de la variance des erreurs est :
\[\hat{\sigma}^2_\varepsilon = \frac{SCR}{n-m-1},\]
la matrice de variance covariance est :
\[\mathbb{V}(\hat{\beta}) = \hat{\sigma}^2_\varepsilon \left( \boldsymbol X^t \boldsymbol X \right)^{-1}\]

> (hatSigma2_u <- SCR / (n-m-1))
[1] 0.6678786
> (varcov <- hatSigma2_u * solve(t(X) %*% X))
                x1            x2            x3            x4          cste
x1    7.622460e-11 -1.256183e-09 -5.149958e-08  3.693662e-09 -3.444875e-07
x2   -1.256183e-09  3.304204e-07 -3.211509e-06 -2.812140e-07  2.823964e-05
x3   -5.149958e-08 -3.211509e-06  9.983711e-05 -1.378010e-06  1.014358e-03
x4    3.693662e-09 -2.812140e-07 -1.378010e-06  7.722665e-07 -4.091131e-04
cste -3.444875e-07  2.823964e-05  1.014358e-03 -4.091131e-04  3.929010e-01
# Les erreurs-types sont les racines des elements diagonaux
> (hatSigmaBetas <- sqrt(diag(varcov)))
          x1           x2           x3           x4         cste 
8.730670e-06 5.748221e-04 9.991852e-03 8.787869e-04 6.268181e-01 

Le test de significativité pour chaque coefficient \(\beta\) est le suivant :
\begin{align*}
\begin{cases}
H_0 : \beta = 0\\
H_1 : \beta \ne 0
\end{cases}.
\end{align*}
Il s'appuie sur la statistique :
\[T = \frac{\beta - 0}{\hat{\sigma}_{\hat{\beta}}} \sim \mathcal{S}t(n-m-1),\]
où \(\hat{\sigma}_{\hat{\beta}}\) est l'estimation de l'écart-type de l'estimateur du paramètre \(\beta\).

La règle de décision est la suivante : si la valeur absolue de la statistique observée est supérieure à la valeur théorique de la Student à \((n-m-1)\) degrés de libertés, pour un risque \(\alpha\) donné, on rejette au seuil de \(\alpha\) l'hypothèse nulle en faveur de l'hypothèse alternative.

Admettons qu'on choisisse (pour être original) un risque de première espèce de \(\alpha=5\%\).

> (tObs <- hatBeta / hatSigmaBetas)
           [,1]
x1    2.3389432
x2   -0.8708666
x3    2.5014857
x4    4.7344612
cste  3.9186711
> # On doit comparer les t_obs a la valeur tabulee d'une
> # Student a (n-m-1) d.d.l. pour un risque de premiere espece
> # alpha = 5%
> (T_tab <- qt(p = 1-0.05/2, df = n-m-1))
[1] 2.055529
> ifelse(abs(tObs) > T_tab, "Rejet H_0", "Non rejet H_0")
     [,1]           
x1   "Rejet H_0"    
x2   "Non rejet H_0"
x3   "Rejet H_0"    
x4   "Rejet H_0"    
cste "Rejet H_0" 

Ainsi, au seuil de \(5\%\), on rejette l'hypothèse de nullité statistique du coefficient associé à chaque coefficient, excepté celui associé à la variable \(x_2\). La p-value (probabilité d'obtenir une valeur au moins aussi grande de la statistique observée, si l'hypothèse nulle est vraie) associée à chaque test est la suivante :

> 2*pt(q = abs(tObs), df = n-m-1, lower.tail=FALSE)
             [,1]
x1   2.729732e-02
x2   3.917972e-01
x3   1.899350e-02
x4   6.774484e-05
cste 5.777629e-04

Ensuite, on peut effectuer le test de globalité de Fisher, qui est le suivant :
\begin{align*}
\begin{cases}
H_0 : \beta_1 = \beta_2 = \beta_3 = \beta_4 = 0\\
H_1 : \textrm{au moins un des \(\beta\) est différent de \(0\)}
\end{cases}
\end{align*}

La statistique de test est la suivante :
\[F = \frac{R^2/m}{(1-R^2)/(n-m-1)} \sim \mathcal{F}(m,n-m-1).\]

À nouveau, on doit comparer la valeur calculée à la valeur théorique. Si la valeur calculée dépasse la valeur théorique, on rejette l'hypothèse nulle, au seuil donnée. Gardons le seuil de \(\alpha=5\%\) :

> # La valeur calculee
> (F_obs <- (R2 / m) / ((1-R2) / (n-m-1)))
[1] 136.5413
> # La valeur tabulee est :
> (F_tab <- qf(p = 1-0.05, df1 = m, df2 = n-m-1))
[1] 2.742594
> # p-value
> (1 - pf(q = F_obs, df1 = m, df2 = n-m-1))
[1] 0

On rejette donc \(H_0\) au seuil de \(5\%\).

Retrouvons à présent ces résultats à l'aide de deux lignes de code R :

> df <- data.frame(y = y, x1 = x1, x2 = x2, x3 = x3, x4 = x4)
> summary(reg <- lm(y~., data = df))

Call:
lm(formula = y ~ ., data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5677 -0.6704  0.1183  0.5283  1.4361 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.456e+00  6.268e-01   3.919 0.000578 ***
x1           2.042e-05  8.731e-06   2.339 0.027297 *  
x2          -5.006e-04  5.748e-04  -0.871 0.391797    
x3           2.499e-02  9.992e-03   2.501 0.018993 *  
x4           4.161e-03  8.788e-04   4.734 6.77e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8172 on 26 degrees of freedom
Multiple R-squared:  0.9546,	Adjusted R-squared:  0.9476 
F-statistic: 136.5 on 4 and 26 DF,  p-value: < 2.2e-16

Dans la fonction lm, le point indique qu'on souhaite régresser \(y\) sur toutes les autres variables de la data.frame. On peut écrire, de manière équivalente :

> summary(reg <- lm(y ~ 1 + x1 + x2 + x3 + x4, data = df))

Faisons comme si le modèle était valide, et donnons une indication de lecture des coefficients. On lit que le coefficient associé à la variable \(x_1\) est \(2.042 \times 10^{-5}\), ce qui signifie que lorsque \(x_1\) diminue d'une unité, \(y\) diminue de \(2.042 \times 10^{-5}\) unités, toutes choses égales par ailleurs. Le coefficient associé à \(x^2\) n'est pas significativement différent de zéro. On ne l'interprète pas.

En fait, on peut voir que \(x_2\) est fortement corrélé aux autres variables explicatives :

> cor(df)
           y        x1        x2        x3        x4
y  1.0000000 0.8911104 0.9409920 0.9526249 0.8638623
x1 0.8911104 1.0000000 0.8977790 0.9351708 0.6349611
x2 0.9409920 0.8977790 1.0000000 0.9625134 0.8378676
x3 0.9526249 0.9351708 0.9625134 1.0000000 0.7798228
x4 0.8638623 0.6349611 0.8378676 0.7798228 1.0000000

On abordera ce problème lors du prochain exercice.

One thought on “[L3 Eco-Gestion] Régression linéaire multiple avec R

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *

Time limit is exhausted. Please reload CAPTCHA.