Avec les L3 Éco-Gestion, on commence à mettre les mains dans le cambouis pour effectuer les premières régressions, avec Excel.
Je souhaite proposer de revenir sur les exercices des TP, en fournissant une alternative à l’utilisation d’Excel, à savoir R. Pourquoi ? J’admets qu’utiliser un tableur puisse aider dans la compréhension des calculs effectués pour faire une régression par moindres carrés ordinaires, mais je ne pense pas que ce soit la seule façon de le faire, et à mon sens, R offre également la possibilité de manier les données, et de voir comment fonctionnent les différents mécanismes de la régression.

Il existe de nombreuses introductions à R. Entre autres, en voici quelques unes qui sont faciles d’accès :

Alors, dans l’exercice, on a seulement 15 observations de deux variables \(x\) et \(y\). On souhaite expliquer les variations de \(y\) par celles de \(x\). Voici les données :

> y <- c(3,4,4,7,9,10,10,12,12,17,14,34,24,21,37)
> y
 [1]  3  4  4  7  9 10 10 12 12 17 14 34 24 21 37
> x <- c(40,46,53,58,65,70,76,82,88,94,101,106,111,117,123)
> x
 [1]  40  46  53  58  65  70  76  82  88  94 101 106 111 117 123

Lorsque la ligne commence par un chevron (>), cela signifie qu’il s’agit d’une instruction passée à R. On peut créer une data.frame (cela correspond à une matrice à double entrée dans laquelle les lignes sont les observations et les colonnes servent à séparer les variables) :

> df <- data.frame(y,x)
> df
    y   x
1   3  40
2   4  46
3   4  53
4   7  58
5   9  65
6  10  70
7  10  76
8  12  82
9  12  88
10 17  94
11 14 101
12 34 106
13 24 111
14 21 117
15 37 123

Puisqu’on veut modéliser la relation entre \(y\) et \(x\) linéairement, regardons si à l’oeil nu cette idée ne paraît pas saugrenue.

> plot(x = x, y = y, t = "p", main = "y en fonction de x")

Alternativement, on peut utiliser (puisque c’est à la mode) le package ggplot2 pour réaliser ce graphique. Lors de la première utilisation, il faut installer la librairie pour pouvoir utiliser les fonctions qui la composent.

install.packages("ggplot2")

Dans tous les autres cas, il faut charger la bibliothèque (une fois par session suffit) :

> library(ggplot2)
> ggplot() + geom_point(data = df, aes(x = x, y = y)) +
+ ggtitle("y en fonction de x")

Cette représentation ne va pas contre l’idée d’une relation linéaire entre \(y\) et \(x\). On peut toutefois se demander ce qui se passe avec les deux points qui ont l’air éloigné. Dans ce contexte, on ne peut rien dire, vu qu’on ne sait ni ce que sont les données, ni comment elles ont été obtenues. Nous allons faire avec.

Le modèle que nous allons estimer est le suivant :
\[y_i = \alpha x_i + \beta + u_i, \quad i = 1,2, \ldots 15\]
avec \(u_i\) un terme d’erreur que l’on suppose normal, d’espérance nulle et de variance \(\sigma^2_u\) identique pour tous les individus.

Nous cherchons les estimateurs \((\hat{\alpha},\hat{\beta})\) par moindres carrés ordinaires (MCO) de \((\alpha, \beta)\). Le cours nous donne :
\begin{align*}
\begin{cases}
\hat{\alpha} = \frac{\sum(x_i y_i) – n \bar{x} \bar{y}}{\sum(x_i^2) – n \bar{x}^2}\\
\hat{\beta} = \bar{y} – \hat{\alpha} \bar{x}
\end{cases}.
\end{align*}

> n <- length(y)
> (alphaChap <- (sum(x*y) - n * mean(x) * mean(y)) / (sum(x^2) - n * mean(x)^2))
[1] 0.3485219
> (betaChap <- mean(y) - alphaChap * mean(x))
[1] -14.04546

On a donc
\begin{align*}
\begin{cases}
\hat{\alpha} = 0.3485219\\
\hat{\beta} = -14.04546
\end{cases}.
\end{align*}

Les valeurs estimées, \(\hat{y}_i\) sont obtenues par :
\[\hat{y}_i = \hat{\alpha} x_i + \hat{\beta}.\]

> (yChap <- alphaChap * x + betaChap)
 [1] -0.1045872  1.9865443  4.4261978  6.1688073  8.6084608 10.3510703 12.4422018 14.5333333 16.6244648
[10] 18.7155963 21.1552497 22.8978593 24.6404689 26.7316004 28.8227319

Il est facile de représenter les valeurs estimées et celles observées sur un graphique :

> ggplot() + geom_line(aes(x = x, y = y, col = "Obs.")) +
+   geom_point(aes(x = x, y = y, col = "Obs.")) +
+   geom_line(aes(x = x, y = yChap, col = "Estim.")) +
+   geom_point(aes(x = x, y = yChap, col = "Estim.")) +
+   ggtitle("Valeurs observees et valeurs estimees par MCO") +
+   scale_colour_discrete("Serie")

Les résidus sont obtenus par \(e_i = y_i - \hat{y}_i\) :

> (residus <- y - yChap)
 [1]  3.1045872  2.0134557 -0.4261978  0.8311927  0.3915392 -0.3510703 -2.4422018 -2.5333333 -4.6244648
[10] -1.7155963 -7.1552497 11.1021407 -0.6404689 -5.7316004  8.1772681

La somme des carrés des résidus, \(SCR = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2\) vaut :

> (SCR <- sum(residus^2))
[1] 326.1369

La somme des carrés totale, \(SCT = \sum_{i=1}^{n} (y - \bar{y})^2\), avec \(\bar{y} = n^{-1}\sum_{i=1}^{n} y_i\), vaut quant à elle :

> (SCT <- sum((y - mean(y))^2))
[1] 1517.733

Ces deux valeurs permettent d'obtenir la somme des carrés expliquée : \(SCE = SCT - SCR\) :

> (SCE <- SCT - SCR)
[1] 1191.596

On peut également calculer le coefficient de détermination (\(R^2\)) à l'aide de la SCR et de la SCT :
\[R^2 = 1 - \frac{SCR}{SCT}\]

> (R2 <- 1 - SCR/SCT)
[1] 0.7851158.

Ainsi, 78.51 pour cent des variations de \(y\) sont expliquées par notre modèle.

La valeur du coefficient de détermination ajusté s'obtient à l'aide de la formule :
\[R^2_a = 1 - \frac{n-1}{n-m-1} \times (1- R^2),\]
avec \(m\) le nombre de variables explicatives, soit une seule dans notre exemple.

> (R2a <- 1 - (n-1)/(n-1-1) * (1-R2))
[1] 0.7685863.

La variance des résidus peut s'obtenir en utilisant la formule :
\[\hat{\sigma}^2_u = \frac{SCR}{n-2}.\]

> (hatSigma2_u <- SCR / (n-2))
[1] 25.08745

ou encore à l'aide de :
\[\hat{\sigma}^2_u = \frac{1-R^2}{n-2} \times \left(\sum_{i=1}^{n}y^2 - n \times \bar{y}^2\right).\]

> (1-R2)/(n-2) * (sum(y^2) - n * mean(y)^2)
[1] 25.08745

La variance estimée des estimateurs de \(\alpha\) et \(\beta\) sont données par :
\begin{align*}
\begin{cases}
\hat{\sigma}^2_{\hat{\alpha}} & = \frac{\hat{\sigma}^2_u}{\sum_{i=1}^{n} x_i^2 - n \bar x^2}\\
\hat{\sigma}^2_{\hat{\beta}} & = \bar x^2 \hat{\sigma}^2_{\hat{\alpha}} + \frac{\hat{\sigma}^2_u}{n}
\end{cases}
\end{align*}

> (hatSigma2_alphaChap <- hatSigma2_u / (sum(x^2) - n * mean(x)^2))
[1] 0.002557335
> (hatSigma2_betaChap <- mean(x)^2 * hatSigma2_alphaChap + hatSigma2_u / n)
[1] 18.86802

La racine carrée de ces estimations de variances permet par exemple, de calculer les statistiques observées de Student, dans le cadre des tests de nullité de coefficients.

> (hatSigma_alphaChap <- sqrt(hatSigma2_alphaChap))
[1] 0.0505701
> (hatSigma_betaChap <- sqrt(hatSigma2_betaChap))
[1] 4.343733

On souhaite effectuer le test suivant :
\begin{align*}
\begin{cases}
H_0 : \alpha = 0\\
H_1 : \alpha \ne 0,
\end{cases}
\end{align*}

La statistique de test est \(T = \frac{\hat{\alpha} - 0}{\hat{\sigma}_{\hat{\alpha}}}\sim \mathcal{S}t(n-2)\). On va devoir comparer la valeur observée \(t_{obs}\) à la valeur théorique d'une loi de Student \(t\). On cherche alors les quantiles \(-t\) et \(t\) tels que \(\mathbb{P}(-t < T < t) = 1-p\).

Ainsi, si la valeur absolue de \(t_{obs}\) dépasse celle de \(t\), on rejette l'hypothèse nulle, au seuil de \(p=5\%\). La valeur observée est :

> (tObs_alpha <- (alphaChap - 0) / hatSigma_alphaChap)
[1] 6.891858

On compare à la valeur théorique :

> qt(p = 1-0.05/2, df = n-2)
[1] 2.160369

On a donc \(\mid t_{obs}\mid > t_{p/2}(15-2)\), on rejette donc l'hypothèse nulle avec un risque de première espèce de \(p=5\%\). Le coefficient \(\alpha\) est, au seuil de \(5\%\), différent de la valeur zéro.

Une alternative pour conclure sur le test, est de s'intéresser à la p-value, c'est à dire la probabilité d'obtenir une statistique de test au moins aussi grande que celle que nous avons observé, si jamais l'hypothèse nulle est vraie. Si cette probabilité d'apparition est inférieure à \(5\%\), on rejette l'hypothèse nulle, au seuil de \(p=5\%\).
La p-value est alors donnée par \(2 \times \mathbb{P}(T > \mid t_{obs} \mid)\) :

> (2 * pt(abs(tObs_alpha), n-2, lower.tail=FALSE))
[1] 1.098635e-05

La probabilité de rejeter à tort l'hypothèse nulle est de \(1.098635\times 10^{-05}\).

Rapidement, pour \(\beta\), on a :

> (tObs_beta <- (betaChap - 0) / hatSigma_betaChap)
[1] -3.233501
> (2 * pt(abs(tObs_beta), 13, lower.tail=FALSE))
[1] 0.00653243

On en tire donc les mêmes conclusions.

Il est possible d'effectuer le test de Fisher de significativité globale, construit à partir de la statistique suivante :
\[F = \frac{R^2 / m}{(1-R^2) / (n-m-1)} \sim \mathcal{F}(m, n-m-1).\]

Cette fois, le test est unilatéral. Si la valeur observée de la statistique est supérieure à la valeur tabulée d'une Fisher à \((m, n-m-1)\) degrés de liberté, on rejette l'hypothèse nulle de nullité simultanée de tous les coefficients des variables explicatives, au seuil donné.

La valeur observée vaut :

> (Fobs <- (R2/1) / ((1-R2) / (n-2)))
[1] 47.4977

que l'on compare à la valeur théorique \(f_{1-p}(1,13)\) :

> qf(p = 1-0.05, df1 = 1, df2 = n-2)
[1] 4.667193

La valeur observée étant largement supérieure à la valeur théorique, on rejette l'hypothèse nulle au seuil de \(p=5\%\). On considère que la régression n'est pas mauvaise (mais on garde à l'esprit que le nombre de données à notre disposition est très faible, et que les tests proposés ici sont valables asymptotiquement).

Comme pour les tests bilatéraux de nullité des coefficients, on peut dans le cadre de ce test unilatéral calculer la p-value :

> 1-pf(q = Fobs, df1 = 1, df2 = n-2)
[1] 1.098635e-05
[code][/code]

La probabilité de rejeter à tort l'hypothèse nulle (rejeter \(H_0\) alors que \(H_0\) est vraie) est quasi nulle ici.

Enfin, on peut calculer des intervalles de confiances pour les paramètres estimés. Ces derniers sont définis comme suit :
\begin{align*}
\widehat{\textrm{I.C.}_{\alpha}(1-p)} & = \left[ \hat{\alpha} \pm t_{\alpha / 2, n-2} \times \hat{\sigma}_{\hat{\alpha}} \right]\\
\widehat{\textrm{I.C.}_{\beta}(1-p)} & = \left[ \hat{\beta} \pm t_{\alpha / 2, n-2} \times \hat{\sigma}_{\hat{\beta}} \right].
\end{align*}

> (alphaChap_inf <- alphaChap - qt(p = 1-0.05/2, df=n-2) * hatSigma_alphaChap)
[1] 0.2392719
> (alphaChap_sup <- alphaChap + qt(p = 1-0.05/2, df=n-2) * hatSigma_alphaChap)
[1] 0.457772
> 
> (betaChap_inf <- betaChap - qt(p = 1-0.05/2, df=n-2) * hatSigma_betaChap)
[1] -23.42953
> (betaChap_sup <- betaChap + qt(p = 1-0.05/2, df=n-2) * hatSigma_betaChap)
[1] -4.661399

Maintenant qu'on a vu comment calculer tous ces différents éléments, on peut utiliser deux lignes de code pour obtenir la même chose, à l'aide de routines déjà établies sous R.

> summary(reg <- lm(y~x))

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.1552 -2.4878 -0.4262  1.4223 11.1021 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -14.04546    4.34373  -3.234  0.00653 ** 
x             0.34852    0.05057   6.892  1.1e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.009 on 13 degrees of freedom
Multiple R-squared:  0.7851,	Adjusted R-squared:  0.7686 
F-statistic:  47.5 on 1 and 13 DF,  p-value: 1.099e-05

> confint(reg, level = 0.95)
                  2.5 %    97.5 %
(Intercept) -23.4295283 -4.661399
x             0.2392719  0.457772

Le lecteur désireux d'extraire les informations contenues dans les résultats de la régression est invité à consulter le code qui suit.

# On peut y lire les estimateurs de \alpha et \beta
# Les estimations de leur variance
# Les t_obs ainsi que les p-values des tests de nullite de chaque coefficient
# Le R^2 et le R^2_a
# Le F_obs ainsi que la p-value du test

# On peut obtenir de nombreuses informations supplementaires :
names(reg)
names(summary(reg))

# Par exemple, pour obtenir l'estimation de la variance des aleas, on peut faire :
summary(reg)$sigma^2

# Ou pour la somme des carres des residus :
sum(residuals(reg)^2)
# Ou de maniere similaire
sum(reg$residuals^2)
sum(summary(reg)$residuals^2)

# Les yChap sont ici :
reg$fitted.values
fitted(reg)

# Le R^2
summary(reg)$r.squared

# Le R^2_a :
summary(reg)$adj.r.squared

# Pour construire les intervalles de confiance a 95%
# pour les estimateurs de \alpha et \beta
confint(reg, level = 0.95)

Sinon, La correction de cet exercice sous Excel est aussi en ligne.

2 thoughts on “[L3 Eco-Gestion] Régression linéaire simple 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.