# 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

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 :

yi=β0+β1xi1+β2xi2+ϵi

où :

• y est la variable à expliquer (ici : sale.price)
• x1 et x2 sont les variables explicatives (ici : la surface et le nombre de chambres),
• ϵi est le terme d’erreur, supposé être indépendants et identiquement distribués selon une loi Normale ϵi∼N(0,σ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 :

{ε∼N(0n,σ2In) ici n=15rg(X)=p=3

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 sansbiais_ et de varianceminimale_, c’est ce que nous assure le théorème de Gauss-Markov. 

Sans biais : Soit ˆγ un estimateur de θ, ˆγ est un estimateur sans biais de θ si E[ˆγ]=θ

Variance minimale : Soit ˆγ un estimateur sans biais de θ, on dit que ˆγ est de variance minimale si pour autre estimateur ˆα sans biais de θ on a : V[ˆγ]≤V[ˆα]

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 S1≤S2 signifie que S=(S2−S1) est une matrice symétrique réelle positive, c’est-à-dire que pour tout vecteur x on a xtS1x≤xtS2x. Cela revient aussi à dire que les valeurs propres de S sont toutes supérieures ou égales à 0.
Ainsi montrer que V[ˆγ]≤V[ˆα] revient à montrer que V[ˆγ]−V[ˆα] est symétrique définie positive c’est-à-dire que pour tout vecteur x on a : xtV[ˆγ]x≤xtV[ˆα]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 : (ˆβ0,ˆβ1)=argminβ0,β1n∑i=1ˆϵ2i=argminβ0,β1n∑i=1(yi−β0−β1x1i−β2x2i)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 : {sale.price}=70.75+0.188×{area}+ˆϵ

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") 

#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 : {sale.price}=−141.76+0.14×{area}+58.32×{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 (H0) : βk=0 pour k∈[[1;3]] vs (H1) : βk≠0 pour k∈[[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 R2=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 R2 défini par R2=1−SSRSST

représente une manière d’évaluer l’adéquation entre le modèle et les données observées. On a R2∈[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 R2adjusted 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=ˆϵiˆσ(i)√1−hii=yi−ˆyiˆσ(i)√1−hii∼Tn−p−1.

Où hii 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 ˆσ(i) est l’estimateur de σ 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.