Rappel :

Le but de la régression linéaire simple est de chercher une fonction \(f\) linéaire telle que \(y_i \approx f(x_i)\). Pour définir \(\approx\), il faut se donner un critère quantifiant la qualité de l’ajustement de la fonction \(f\) aux données. Ainsi, une étude de régression simple débute généralement par un tracé des observations \((x_i, y_i)\), \(i = 1, ..., n\). Où les \(y_i\) représentent les données de la variable à expliquer et les \(x_i\) représentent les données de la variable explicative. Cette première représentation permet de savoir si le modèle linéaire est pertinent. Le graphique suivant représente trois nuages de points différents :

# Générer des données aléatoires
set.seed(123)
n <- 100
x <- seq(1, 10, length.out = n)

# Définir la disposition des graphiques
par(mfrow = c(1, 3))

# Premier graphique : nuage de points avec une forme sinusoïdale
y_sinusoide <- sin(x) + rnorm(n, mean = 0, sd = 0.5)
plot(x, y_sinusoide, main = "Graphique 1", xlab = "X", ylab = "Y")

# Deuxième graphique : nuage de points avec une forme sigmoïdale
y_sigmoide <- 1 / (1 + exp(-(x - 5))) + rnorm(n, mean = 0, sd = 0.3)
plot(x, y_sigmoide, main = "Graphique 2", xlab = "X", ylab = "Y")

# Troisième graphique : nuage de points propice à la régression linéaire
y_lineaire <- 2 * x + 3 + rnorm(n, mean = 0, sd = 1)
plot(x, y_lineaire, main = "Graphique 3", xlab = "X", ylab = "Y")

Nous avons ici simuler des données fictives avec une forme sinusoïdale et sigmoïdale. Au vu des graphiques, il semble inadéquat de proposer une régression linéaire pour les 2 premiers graphiques, le tracé présentant une forme sinusoïdale ou sigmoïdale. Par contre, la modélisation par une droite de la relation entre \(X\) et \(Y\) pour le dernier graphique semble correspondre à une bonne approximation de la liaison.

A noter que dans le cas du graphique 1 et 2, d’autres techniques de modélisations existent comme des méthodes de régression linéaires de type polynomiale, exponentielle, logarithmique…

Régression linéaire multipe avec R :

Pour illustrer, on utilise le jeu de données « housingprices » (issu du package DAAG), composé de quinze observations et trois variables qui sont :

  1. sale.price := prix de vente de la maison (en milliers de dollards australiens)
  2. area = surface au sol de la maison
  3. bedrooms = nombre de chambres dans la maison
# Charge le package DAAG

library(DAAG)

# Charge les données 
data(houseprices)

houseprices
##    area bedrooms sale.price
## 9   694        4      192.0
## 10  905        4      215.0
## 11  802        4      215.0
## 12 1366        4      274.0
## 13  716        4      112.7
## 14  963        4      185.0
## 15  821        4      212.0
## 16  714        4      220.0
## 17 1018        4      276.0
## 18  887        4      260.0
## 19  790        4      221.5
## 20  696        5      255.0
## 21  771        5      260.0
## 22 1006        5      293.0
## 23 1191        6      375.0

Ce dataframe contient 15 données des informations sur les prix des maisons ainsi que diverses caractéristiques associées à ces maisons.

Objectif :

On peut vouloir expliquer le prix de vente des maisons, en fonction de leur surface en présumant que plus la surface est élevée, et plus le prix de vente sera élevé. Il s’agit là d’une regression dite “simple” car elle ne comporte qu’une seule variable explicative. De plus, on peut aussi supposer que le nombre de chambres influence le prix de la maison à la hausse ; il s’agira là d’une regression “multiple” avec deux variables explicatives.

L’objectif est d’évaluer si chacune des deux variables influence le prix, et, si tel est le cas, de tenter de quantifier cet effet.

Représentation graphique :

# Nuage des 15 points entre 'area' et 'sale.price'
plot(houseprices$area, houseprices$sale.price, 
     xlab = 'Area', ylab = 'Sale Price', 
     main = 'Nuage de points Area vs Sale Price')

La régression linéaire assure que la relation entre les variables explicatives et la variable à expliquer (variable numérique continue) va être linéaire, du type :

\(y_{i} = \beta_0 + \beta_1 x_{i_1} + \beta_2 x_{i_2} + \epsilon_{i}\)

où :

  • \(y\) est la variable à expliquer (ici : sale.price)
  • \(x_1\) et \(x_2\) sont les variables explicatives (ici : la surface et le nombre de chambres),
  • \(\epsilon_i\) est le terme d’erreur, supposé être indépendants et identiquement distribués selon une loi Normale \(\epsilon_i \sim \mathcal{N}(0,\sigma^2)\)

Mieux que les expressions des estimateurs et celles de leurs variances, on aimerait connaître leurs lois : ceci permettrait par exemple d’obtenir des régions de confiance et d’effectuer des tests d’hypothèses. Dans cette optique, il faut bien entendu faire une hypothèse plus forte sur notre modèle, à savoir préciser la loi des erreurs. Les hypothèses sont alors :

\[ \left\{ \begin{array}{c} \varepsilon \sim N(0_n, \sigma^2I_{n})\ ici\ n=15 \\ rg(X)=p=3\\ \end{array} \right. \]

Ainsi la méthode des moindres carrés ordinaires (OLS pour Ordinary Least Square) permet d’obtenir des estimateurs BLUE (Best Linear Unbiased Estimators), c’est-à-dire qui sont \(\underline{sans \; biais}\) et de \(\underline{variance \; minimale}\), c’est ce que nous assure le théorème de Gauss-Markov. 

Sans biais : Soit \(\hat\gamma\) un estimateur de \(\theta\), \(\hat\gamma\) est un estimateur sans biais de \(\theta\) si \(\mathbb{E}[\hat\gamma] = \theta\)

Variance minimale : Soit \(\hat\gamma\) un estimateur sans biais de \(\theta\), on dit que \(\hat\gamma\) est de variance minimale si pour autre estimateur \(\hat\alpha\) sans biais de \(\theta\) on a : \(\mathbb{V}[\hat\gamma] \leq \mathbb{V}[\hat\alpha]\)

Dans le cas de la régression linéaire multiple on peut en fait définir une relation d’ordre partielle entre matrices symétriques réelles : dire que \(S_1 \leq S_2\) signifie que \(S= (S_2 - S_1)\) est une matrice symétrique réelle positive, c’est-à-dire que pour tout vecteur \(x\) on a \(x^{t}S_{1}x \leq x^{t}S_{2}x\). Cela revient aussi à dire que les valeurs propres de S sont toutes supérieures ou égales à 0.
Ainsi montrer que \(\mathbb{V}[\hat\gamma] \leq \mathbb{V}[\hat\alpha]\) revient à montrer que \(\mathbb{V}[\hat\gamma]-\mathbb{V}[\hat\alpha]\) est symétrique définie positive c’est-à-dire que pour tout vecteur \(x\) on a : \(x^{t}\mathbb{V}[\hat\gamma]x \leq x^{t}\mathbb{V}[\hat\alpha]x\) .

Les coefficients estimés sont issus de la minimisation de la somme des carrés des résidus entre les valeurs prédites et les valeurs observées, formellement : \[(\hat{\beta}_0, \hat{\beta}_1) = \underset{\beta_0, \beta_1}{\arg\min} \sum_{i=1}^{n} \hat{\epsilon}_i^2 = \underset{\beta_0, \beta_1}{\arg\min} \sum_{i=1}^{n} (y_i - \beta_0 - \beta_1 x1_i - {\beta}_2 x2_i)^2\]

Régression linéaire simple

# y = b0 + b1*x1
# Variable à expliquer : y = prix de vente de la maison
# Une variable explicative : x1 = Surface au sol de la maison

pricereg<-lm(sale.price ~ area, data=houseprices)
summary(pricereg)
## 
## Call:
## lm(formula = sale.price ~ area, data = houseprices)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -92.499 -19.302   2.406  28.019  80.607 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  70.7504    60.3477   1.172   0.2621  
## area          0.1878     0.0664   2.828   0.0142 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 48.18 on 13 degrees of freedom
## Multiple R-squared:  0.3809, Adjusted R-squared:  0.3333 
## F-statistic: 7.997 on 1 and 13 DF,  p-value: 0.01425

On a donc l’équation de la droite de régression : \[ \text{{sale.price}} = 70.75 + 0.188 \times \text{{area}} + \hat{\epsilon} \]

que l’on peut tracer avec le nuage de points :

# plot : "vraies" valeurs et droite de regression
plot(sale.price ~ area, data=houseprices)
abline(pricereg, col = "red")

Régression linéaire multiple :

On considère deux (ou plus) variables explicatives. On souhaite donc estimer les coefficients du modèle : \[ \text{{sale.price}} = \beta_0 + \beta_1 \times \text{{area}} + \beta_2 \times \text{{bedroom}} + \epsilon \]

#variables explicatives : x1 = Surface au sol de la maison, x2 = nombre de chambres
pricereg2<-lm(sale.price ~ area + bedrooms , data=houseprices)
summary(pricereg2)
## 
## Call:
## lm(formula = sale.price ~ area + bedrooms, data = houseprices)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -80.897  -4.247   1.539  13.249  42.027 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -141.76132   67.87204  -2.089  0.05872 . 
## area           0.14255    0.04697   3.035  0.01038 * 
## bedrooms      58.32375   14.75962   3.952  0.00192 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.06 on 12 degrees of freedom
## Multiple R-squared:  0.731,  Adjusted R-squared:  0.6861 
## F-statistic:  16.3 on 2 and 12 DF,  p-value: 0.0003792

Interprétation :

La droite de régression est donnée par l’équation : \[ \text{{sale.price}} = -141.76 + 0.14 \times \text{{area}} + 58.32 \times \text{{bedrooms}} \]

Les coefficients associés aux variables “area” et “bedrooms” sont significatifs (respectivement à 95% et à 99%, cf. sortie précédente) d’après le test de Student qui a pour hypothèses \((H_0)\) : \(\beta_k = 0\ pour\ k\in [\![1;3]\!]\) vs \((H_1)\) : \(\beta_k\ne 0\ pour\ k\in [\![1;3]\!]\) , en effet on voit que dans Pr(>|t|) ces valeurs sont inférieures à 0.05, ces coefficients sont donc significatifs, on peut donc les interpréter. Toutes choses égales par ailleurs, une chambre supplémentaire augmente par exemple le prix de la maison d’environ 58 mille dollars, et 100 unités de surface (inconnue) supplémentaires vont l’augmenter de 14 mille dollars, toutes choses égales par ailleurs. De plus on voit que le \(R^2=0,731\) notre modèle de régression explique donc 73,1% de de la variance des données.

Commandes utiles :

La commande summary() affichera des informations telles que les coefficients estimés, les erreurs standards, les valeurs t et p associées, le coefficient de détermination (R-squared), et d’autres statistiques pertinentes pour évaluer la qualité du modèle.

Par exemple le \(R^2\) défini par \[ R^2 = 1 - \frac{SSR}{SST} \] représente une manière d’évaluer l’adéquation entre le modèle et les données observées. On a \[ R^2 \in [0, 1] \] et plus sa valeur est proche de 1 plus l’adéquation sera forte. Cependant la valeur est influencée par le nombre de variables explicatives incluses dans la régression. Ainsi le \(R^{2}_{adjusted}\) va tenir compte de ce nombre est sera plus correct. La statistique de Fisher est une mesure de l’adéquation globale du modèle. Elle évalue si le modèle de regression dans son ensemble est significativement différent d’un modèle sans aucune variable explicative, c’est-à-dire un modèle qui prédit simplement la moyenne de la variable dépendante pour chaque observation. Il est important de choisir un modèle avec statistique F grande et une p-value qui sera faible.

summary(pricereg2)
## 
## Call:
## lm(formula = sale.price ~ area + bedrooms, data = houseprices)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -80.897  -4.247   1.539  13.249  42.027 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -141.76132   67.87204  -2.089  0.05872 . 
## area           0.14255    0.04697   3.035  0.01038 * 
## bedrooms      58.32375   14.75962   3.952  0.00192 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.06 on 12 degrees of freedom
## Multiple R-squared:  0.731,  Adjusted R-squared:  0.6861 
## F-statistic:  16.3 on 2 and 12 DF,  p-value: 0.0003792

La commande coef() permet de n’extraire que les coefficients estimés.

coef(pricereg2)
##  (Intercept)         area     bedrooms 
## -141.7613221    0.1425469   58.3237508

La commande confint()permet d’afficher l’intervalle de confiance à 95% pour les coefficients estimés.

confint(pricereg2)
##                     2.5 %     97.5 %
## (Intercept) -289.64179643  6.1191523
## area           0.04019939  0.2448945
## bedrooms      26.16530118 90.4822003

La commande fitted() permet d’afficher les valeurs prédites :

fitted(pricereg2)
##        9       10       11       12       13       14       15       16 
## 190.4612 220.5387 205.8563 286.2528 193.5973 228.8064 208.5647 193.3122 
##       17       18       19       20       21       22       23 
## 236.6465 217.9728 204.1458 249.0701 259.7611 293.2596 377.9546

La commande predict() permet de prédire la valeur de y (i.e du prix) pour de nouvelles données (des variables explicatives). Seul les deux premiers arguments sont requis : se.fit permet d’afficher l’écart-type de la valeur prédite, et interval et level permettent afficher ici les valeurs de l’intervalle de confiance fixé à 95%. interval = "confidence" pour spécifier que nous voulons des intervalles de confiance pour l’estimation de nos paramètres, et level = 0.95 pour spécifier un niveau de confiance de 95% pour ces intervalles de confiance. Tandis que interval="prediction" c’est pour obtenir des intervalles de prédiction pour les valeurs y pour de nouvelles données.

predict(pricereg, newdata=data.frame(area=800,bedrooms=2), se.fit=TRUE, interval = "prediction", level = 0.95)
## $fit
##        fit      lwr      upr
## 1 220.9719 112.7066 329.2373
## 
## $se.fit
## [1] 13.78228
## 
## $df
## [1] 13
## 
## $residual.scale
## [1] 48.18185

Ainsi, le prix estimé d’une maison dont la surface au sol est de 800 unités serait de 88.92 mille dollars avec un intervalle de prédiction à 95% qui serait \([112.7066 ; 329.2373]\)

La commande resid() permet d’extraire les résidus (Valeur réelle - valeur prédite)

resid(pricereg2)
##           9          10          11          12          13          14 
##   1.5387504  -5.5386516   9.1436820 -12.2527859 -80.8972821 -43.8063735 
##          15          16          17          18          19          20 
##   3.4352904  26.6878118  39.3535454  42.0271931  17.3542452   5.9299058 
##          21          22          23 
##   0.2388861  -0.2596422  -2.9545748

Plusieurs hypothèses sous-jacentes au modèle de regression linéaire portent sur les résidus de la regression : indépendance, homoscédasticité (même variance) et normalité. On peut représenter graphiquement les résidus, afin d’observer si ces conditions sont (plus ou moins) respectées

Exemple d’interprétation des résidus :

# Calcul des résidus
res <- resid(pricereg2)

# Calcul des résidus standardisés
std_res <- res / sd(res)

# Tracé des résidus standardisés
plot(std_res, main = "Résidus Standardisés")


abline(h=0,col="red")

qqnorm(rstudent(pricereg2)); qqline(rstudent(pricereg2), col = "lightblue", lwd = 2)

Le graphique QQ-plot permet de vérifier si les résidus suivent approximativement une distribution normale. Si la distribution des résidus est proche de la droite bleu, cela indique que l’assomption de normalité des résidus est raisonnable. Ici, les queues de distribution semblent plus lourdes qu’une simple gaussienne. Mais ce n’est pas forcément grave : les résidus studentisés sont de type Student, et non gaussiens, et la loi de Student a des queues plus lourdes que la gaussienne, surtout quand le nombre d’observations est petit (ce qui est le cas ici).

De plus la commande rstudent représente les résidus studentisés ils sont utilisés pour évaluer l’influence des observations individuelles sur le modèle de régression. Ils sont également utilisés pour détecter les valeurs aberrantes dans l’ensemble de données. On sait de plus que :

Si la matrice \(X\) est de plein rang et si la suppression de la ligne \(i\) ne modifie pas le rang de la matrice, alors les résidus studentisés par validation croisée vérifie :

\(t^*_i = \frac{\hat{\epsilon}_i}{\hat{\sigma}_{(i)} \sqrt{1 - h_{ii}}} = \frac{y_i - \hat{y}_i}{\hat{\sigma}_{(i)} \sqrt{1 - h_{ii}}}\sim \mathscr{T}_{n - p - 1}.\)

Où \(h_{ii}\) est l’élément de la diagonale de la matrice de projection H qui est la matrice de projection sur l’espace engendré pas les colonnes de X, il quantifie l’effet de l’observation i sur sa propre prédiction et \(\hat{\sigma}_{(i)}\) est l’estimateur de \(\sigma\) dans le modèle linéaire privé de l’observation i.

Autre exemple :

## residus vs. valeurs prédites

fit <- fitted(pricereg2)
plot(fit,res,main="Residuals vs. fitted")
abline(h=0,col="red")

Ici, on remarque que les points sont répartis aléatoirement autour de l’axe horizontal \(y=0\) et qu’il n’y a pas de tendance. Les résidus de ce modèle de régression linéaire sont distribués de manière aléatoire et homogène autour de zéro cela suggère que le modèle de régression linéaire capture correctement la structure des données et que les erreurs du modèle sont essentiellement aléatoires.

Test de significativité globale du modèle :

Sous les hypothèses suivantes :

\(y = X\beta + \epsilon\)

Avec toujours :

\[ \left\{ \begin{array}{c} \varepsilon \sim N(0_{n}, I_{n}\sigma^2) \\ rg(X) = p \end{array} \right. \]

En posant \((H_0)\) : \(\beta_2 = ... = \beta_p = 0\) vs \((H_1)\) : \(\exists i \in [\![2;n]\!]\) tel que \(\beta_i \neq 0\)

La statistique de test associé est :

\(F = \frac{\| \hat{Y} - \bar{y}\mathbf{1} \|^2 / (p - 1)}{\| Y - \hat{Y} \|^2 / (n - p)} \sim F_{n-p}^{p-1}\)

Cela fonctionne si la première colonne de X est l’intercept. C’est le test par défaut renvoyé par R.

Avec \(\mathbb{P}[\text{Rejeter}|H_0 \text{ est vraie}] \leq \alpha\) \(\Leftrightarrow\) \(R_{\alpha} = \{f \in \mathbb{R}_{+} | f \leq f^{p-1}_{n-p}(1-\alpha)\}\) Avec pour p-value : \(p= \mathbb{P}_{H_0}[F>F^{\text{obs}}]\)

Ainsi :

Ce test (F-test) est basé sur la statistique de Fisher présentée en bas de la sortie R. Ici, comme la p-value associée est inférieure à 1%, on peut dire que l’on rejète fortement H0, à savoir le modèle est bien globalement significatif. Si l’on n’avait pas rejeté H0, on n’aurait pas rejeté un modèle où tous les coefficients sont nuls (hors inetercept).

lm(formula = sale.price ~ area + bedrooms, data = houseprices)
## 
## Call:
## lm(formula = sale.price ~ area + bedrooms, data = houseprices)
## 
## Coefficients:
## (Intercept)         area     bedrooms  
##   -141.7613       0.1425      58.3238
summary(lm(formula = sale.price ~ area + bedrooms, data = houseprices))
## 
## Call:
## lm(formula = sale.price ~ area + bedrooms, data = houseprices)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -80.897  -4.247   1.539  13.249  42.027 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -141.76132   67.87204  -2.089  0.05872 . 
## area           0.14255    0.04697   3.035  0.01038 * 
## bedrooms      58.32375   14.75962   3.952  0.00192 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.06 on 12 degrees of freedom
## Multiple R-squared:  0.731,  Adjusted R-squared:  0.6861 
## F-statistic:  16.3 on 2 and 12 DF,  p-value: 0.0003792