R´egression lin´eaire avec R
12 mars 2014
1
Introduction
Nous allons travailler sur le jeu de donn´ees bodyfat issu du package mboost. Pour ce faire, nous entrons les lignes
de commandes suivantes :
> library(mboost)
> data(bodyfat)
Attention : si la commande library retourne une erreur stipulant : “Erreur dans library(”nom”) : aucun package
nomm´e ‘nom’ n’est trouv´e” c’est que le package en question n’est pas encore installer. Ceci peut eˆ tre fait en utilisant la
fonction install.packages("nom du package") (ne pas oublier les guillemets dans cette commande !).
Notez bien que la commande data permet simplement de charger le jeu de donn´ees dans la m´emoire de R. En r´ealit´e,
de la mˆeme mani`ere que l’ouverture de R n’entraˆıne pas le chargement de tous les packages disponibles, le chargement d’un
package par la commande library n’entraˆıne pas le chargement de tous les jeux de donn´ees contenus dans le package.
Face a` un nouveau jeu de donn´ees, il est utile d’obtenir quelques informations e´ l´ementaires telles que : le nombre
d’individus N , le nombre de variables P , la nature de chacune des variables... Les fonctions nrow, ncol et summary
permettent de r´epondre a` ces questions :
> N<-nrow(bodyfat)
> print(N) #Nombre d'individus
[1] 71
> P<-ncol(bodyfat)
> print(P) #Nombre de variables
[1] 10
> summary(bodyfat)
age
Min.
:19.00
1st Qu.:42.00
Median :56.00
Mean
:50.86
3rd Qu.:62.00
Max.
:67.00
kneebreadth
Min.
: 7.200
1st Qu.: 8.600
Median : 9.200
DEXfat
Min.
:11.21
1st Qu.:22.32
Median :29.63
Mean
:30.78
3rd Qu.:39.33
Max.
:62.02
anthro3a
Min.
:2.400
1st Qu.:3.540
Median :3.970
waistcirc
Min.
: 65.00
1st Qu.: 78.50
Median : 85.00
Mean
: 87.38
3rd Qu.: 99.75
Max.
:117.00
anthro3b
Min.
:2.580
1st Qu.:4.060
Median :4.390
1
hipcirc
Min.
: 88.00
1st Qu.: 96.75
Median :103.00
Mean
:105.28
3rd Qu.:111.15
Max.
:132.00
anthro3c
Min.
:2.050
1st Qu.:3.480
Median :3.990
elbowbreadth
Min.
:5.200
1st Qu.:6.200
Median :6.500
Mean
:6.508
3rd Qu.:6.900
Max.
:7.400
anthro4
Min.
:3.180
1st Qu.:5.040
Median :5.530
Mean
: 9.301
3rd Qu.: 9.800
Max.
:11.800
Mean
:3.869
3rd Qu.:4.155
Max.
:4.680
Mean
:4.291
3rd Qu.:4.660
Max.
:5.010
Mean
:3.886
3rd Qu.:4.345
Max.
:4.620
Mean
:5.398
3rd Qu.:5.840
Max.
:6.370
Dans le cadre de cette introduction aux mod`eles lin´eaires avec R, nous ne retiendrons que les 6 premi`eres variables de
ce jeu de donn´ees :
> donnees<-bodyfat[,1:6]
> rownames(donnees)<-paste("ind",1:71)
Notons que la s´election d’une sous-s´election de variables dans un data frame peut aussi se faire en s´electionnant les
variables par leur nom. Ainsi, bodyfat[,c("age","DEXfat")] permet de s´electionner les deux premi`eres variables.
Afin de visualiser graphiquement le jeu de donn´ees, nous utiliserons la fonction ggpairs du package GGally. La
figure obtenue (Figure 1) repr´esente un diagramme de double projection (biplot graphic, en anglais) pour chaque couple de
variable. Des informations suppl´ementaires, comme la distribution de chaque variable et la distribution bivari´ee de chaque
couple de variables peuvent eˆ tre rajout´ees sur le graphique. Notons que la fonction basique de R, pairs permet d’obtenir
un r´esultat e´ quivalent.
> library(GGally)
> ggpairs(donnees, diag=list(continuous="density"),
upper=list(continuous="density"),
lower=list(continuous="smooth"),
axisLabels='show')
` ce stade, plusieurs remarques peuvent eˆ tre formul´ees :
A
— Le jeu de donn´ees contient majoritairement des individus plutˆot aˆ g´es. En effet, la m´ediane est de 56 ans (voir r´esultat
de fonction summary) et la distribution des aˆ ges est tr`es asym´etrique avec un mode au alentour de 60 ans.
— Certaines variables de ce jeu de donn´ees semblent pr´esenter une corr´elation lin´eaire tr`es forte (par exemple la variable DEXfat et la variable waistcirc).
Dans la prochaine partie, nous allons e´ tudier l’influence du tour de taille (variable waistcirc) sur la mesure de la
graisse corporelle (variable DEXfat). Ce choix est naturel au vue de la remarque pr´ec´edente.
2
R´egression lin´eaire simple
Tout d’abord, formalisons notre probl`eme dans un cadre statistique. Nous d´efinissons les variables suivantes :
— soit y = (y1 , ..., yN )0 le vecteur de longueur N = 71 qui contient les valeurs de la variable DEXfat pour chacun
des individus : c’est la variable r´eponse, e´ galement nomm´ee variable endog`ene ou variable a` expliquer,
— soit x1 = (x11 , ..., x1N )0 le vecteur de longueur N = 71 qui contient les valeurs de la variable waistcirc pour
chacun des individus : c’est la variable explicative ou variable exog`ene,
— soient β0 et β1 deux param`etres r´eels inconnus.
Dans l’id´eal, nous voudrions trouver une relation strictement lin´eaire de la forme :
y = β0 + β1 x 1 .
(1)
Si la relation d´ecrite par l’´equation (1) e´ tait v´erifi´ee, elle impliquerait que le diagramme de double projection entre les
variables r´eponse et explicative consiste en une droite parfaite. Nous repr´esentons ces deux variables sur un tel graphique
(Figure 2).
2
age
DEXfat
60
50
40
30
20
10
waistcirc
110
100
90
80
70
kneebreadthelbowbreadth hipcirc
60
50
40
30
●
130
120
110
100
90
7.5
7.0
6.5
6.0
5.5
12
11
10
9
8
7
●
● ●●
●● ●● ●●
●
●
●
●
●
●●●
●
●●
●
●●●●
●
●
●
● ● ●●●●●●
●
●
●●
●
●
●
●
●
●●●● ● ● ●
●
●
●
● ●●●● ●
●
● ●●
●
●
● ● ●●
●●
●
●
●●
●
●
●
●●
●
●
●
● ●
●
●●
● ●
●
● ● ●●
●
●
● ●● ●●
●
●
●
●
● ●●
●● ● ● ●
●
●●
●●
● ●●
●● ●
●
●
●
●●
●
●
●●
●●●
●
●
●
●
●
●
●
●●
●●●
●●
● ●
●● ●
●●●
●
●
●
●●
●
●
●●
●
●●●●●
●
●
●●
●●
●●
●
●
●
●●
●●
● ●
●
●●
●●●
●●
●● ●●
●
●
●●
●● ●
●
●●
●
●
●
●
●●
●●● ●● ●
●
●●
●
●
●
●
●
●
●●
●
●● ●●
●
●
●
●
●
●
●
●
●●●
●●
●●●
●●
●●
●
●
●
●
●
●
●●
●
●●
●
● ●●
●●●
●
●
●
●●●
●
●
●
●●
●
●
●
●
●
●
●
●
● ●●
●●
●
●
● ●
● ●
● ● ● ●●
●●
●
●
●●
●
●●
●●
●
●
●● ●●●
●
●
●
●
●
●
●●
●
●
●
●
●●
● ● ● ●●
●
●
● ● ●
●●
●●
●
●●
●
●●
● ●
●
●●
●
●
●
●
●
●● ●●
●
●
● ●●●●
●
●●
●
●●
●●
●●● ●●
●
●●
●
●
●
●●
●●● ●●●
●●
●
●●
●●
●●
●●●
●● ●●
●
●
●
●●
● ●
●
●
●
●
●
●
●
●
● ●●
●●
●
●
●●
●●
●
●
● ●●●● ●
● ● ●●●
●
●
●
●
●
●
●
●
●
●●
● ●●
● ●
●
●
●
●● ●
●
●●●
●
● ●
●
●
●
●
●
●
●●
●
●
● ●
● ● ●
● ●● ●●
●
●
● ● ●●
●●
● ●
● ●
●● ● ●
●
●
●
●●●●
●●
●
●●●●
●
● ●●
●
●●●
●
●
●●●●●
●
●
●
●●● ●
●
●
●
●●●
● ● ●●
●
● ●●●●
●●
●●●
●●
●
●
●●● ●●
● ●● ●
● ●●
●
●●●
●●
●
●
●
●
●●
●●● ●● ●
●●●
●
● ● ●
●●●
●
●
●
●●
● ●
●
●●
●●
●
●●●
●●
●
●●
●
●●●●
●
●
●
●
●● ●
● ●●
●●●
●
●
● ●● ●
● ●●●●●● ●●
● ●●●
●
●
● ● ●●● ●●
●● ●
●
●● ●
●
● ●
● ● ●●
● ●●●● ●
●
●●
●●● ●
● ●
●
●●
●●
●●●
●● ●●
●●●●●
●
●
●
●
●●●● ●●
●
●
●
●
●
●
●
● ●
●
●
●
●
● ●●● ●●●
●●●
●●●●●
●
●
●
●
●●●
●
●●
●
●●● ●●
●●
●
●
●●●
●
●●
●●
● ●
●●● ● ● ●●
●
●
●
●
●●
●
●
●●●
● ●●
●
●
● ●●
●
● ●
●
●●●
●
●● ●
●● ●
●
●●●
● ●
●● ● ●
●●●●
●●●
●●●● ●
●●●
●
●
●
●●
●
● ●●●●●
●●●●●●
●●
●
●● ●●●●
●●●●● ● ●
●● ●●
● ●●
●●●●●●
●●● ●●
● ●●●●
●
2030405060 102030405060 708090100
110 90100110120130 5.56.06.57.07.57 8 9 101112
age
DEXfat
waistcirc
hipcirc elbowbreadthkneebreadth
F IGURE 1 – Pr´esentation graphique du jeu de donn´ees
> qplot(donnees$waistcirc,donnees$DEXfat,
xlab="Tour de taille",
ylab="Mesure de graisse corporelle",
geom = c("point", "smooth","line"))
Notez que la fonction qplot permet de rajouter sur le graphique des informations suppl´ementaires comme une courbe
liss´ee. Nous constatons que les points semblent s’allonger grossi`erement autour d’une droite. En revanche, il n’existe aucune
droite (ou, autrement dit, aucun couple (β0 , β)) telle que tous les points appartiennent a` cette droite. Le mod`ele propos´e
dans l’´equation (1) n’est donc pas le bon. Pour corriger le mod`ele pr´ec´edent, nous allons eˆ tre oblig´e de consid´erer que les
donn´ees que nous observons contiennent un bruit, une erreur d’ajustement. Notons donc ε = (ε1 , ..., εN )0 le vecteur de
longueur N = 71 qui contient les erreurs al´eatoires de moyenne nulle (ou e´ carts a` la droite de r´egression), et posons le
mod`ele classique de r´egression lin´eaire a` une variable :
y = β0 + β1 x1 + ε.
3
(2)
●
●
Mesure de graisse corporelle
60
●
●
●
●
●● ●
●
●●
40
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
● ●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
20
●
●
●
●
●
●
●
●
●
●
●
70
80
90
100
110
Tour de taille
F IGURE 2 – Mesure de masse corporelle en fonction du tour de taille
2.1
Hypoth`eses relatives au mod`ele lin´eaire
Plusieurs hypoth`eses sont n´ecessaires afin d’obtenir de bonnes propri´et´es des estimateurs (βˆ0 , βˆ1 ) des param`etres du
mod`ele (estimateurs des param`etres sans biais et de variance minimale, par exemple). Nous distinguons deux types d’hypoth`eses : celles que nous appellerons explicites portent sur les termes d’erreur, et celles que nous appellerons implicites
sont directement li´ees a` l’´equation du mod`ele (´equation (??).
Les hypoth`eses explicites, c’est-`a-dire celles qui portent sur les termes d’erreur εi , sont au nombre de trois :
— Ils sont de variance constante, not´ee σ 2 , c’est l’hypoth`ese d’homosc´edasticit´e : en particulier, les erreurs ne sont
pas d´ependantes de la variable exog`ene x1 .
— Ils suivent une loi normale de moyenne nulle et de variance σ 2 .
— Ils sont ind´ependants.
` ces hypoth`eses sur le terme d’erreur, nous devons ajouter deux hypoth`eses suppl´ementaires, qui sont implicites au
A
mod`ele :
— La relation entre la variable exog`ene x1 et la variable endog`ene y est lin´eaire. En particulier, il faudra se m´efier des
relations non-lin´eaires (relation quadratique, exponentielle ou logarithmique par exemple) et des ruptures de pentes.
Nous pr´eciserons ce dernier point dans le paragraphe concern´e.
— La variable exog`ene x1 est mesur´ee sans erreur. En effet, le mod`ele d´ecrit dans l’´equation (2) suppose que la variable
x1 est fixe. Dans la conception originelle du mod`ele lin´eaire, la variable y est e´ galement suppos´ee eˆ tre mesur´ee sans
erreur car le terme al´eatoire ne repr´esente que le d´efaut d’ajustement.
La premi`ere e´ tape de notre travail consiste maintenat a` tester ces hypoth`eses.
2.1.1
V´erification de la condition de lin´earit´e
Le test de corr´elation lin´eaire
4
La premi`ere chose a` faire ici est de tracer la variable r´eponse en fonction de la variable explicative, si cela n’a pas encore
e´ t´e fait (Figure 2). La corr´elation lin´eaire peut eˆ tre calcul´ee grˆace a` la fonction cor :
> cor(donnees$waistcirc,donnees$DEXfat)
[1] 0.8986535
Nous rappelons que la corr´elation lin´eaire varie entre -1 et 1. Lorsque la corr´elation vaut 1, elle indique une corr´elation
positive parfaite entre les donn´ees ; les donn´ees sont parfaitement align´ees le long d’une droite dont le cœfficient directeur
est positif. Lorsqu’elle vaut -1, elle indique une corr´elation n´egative parfaite entre les donn´ees ; ces derni`eres sont alors
parfaitement align´ees le long d’une droite dont le cœfficient directeur est n´egatif.
Attention : si les variables dont nous calculons la corr´elation sont ind´ependantes les unes par apport aux autres, alors le
cœfficient de corr´elation lin´eaire vaudra 0. En revanche, la r´eciproque n’est pas vraie. Regardons l’exemple suivant :
>
>
>
>
#cette commande vous permettra d'obtenir exactement les mˆ
emes r´
esultats
set.seed(1)
exemple<-rnorm(1000,0,1)
cor(exemple,exempleˆ2)
[1] -0.02948134
La corr´elation quadratique entre le vecteur exemple et son carr´e est parfaite ; cependant, la corr´elation lin´eaire entre
les deux vecteurs est proche de z´ero. Dans la Figure 3, nous avons repr´esent´e le vecteur exemple2 en fonction du vecteur
exemple et nous avons ajout´e la droite de r´egression lin´eaire. Le d´efaut de lin´earit´e peut conduire a` une interpr´etation
erronn´ee des r´esultats obtenus. Dans le cas de cet exemple, la Figure 3 montre clairement le d´efaut de lin´earit´e. Voil`a la
commande pour tracer la Figure 3 :
> qplot(exemple,exempleˆ2,xlab="exemple",ylab="exempleˆ2",geom=("point"))+
geom_smooth(method='lm')
Il est possible de tester si la corr´elation observ´ee est significativement diff´erente de 0 ou non. Plus pr´ecis´ement, voici les
hypoth`eses du test :
H0 : la corr´elation lin´eaire entre les deux variables est nulle,
contre
H1 : la corr´elation lin´eaire entre les deux variables est diff´erente de z´ero.
Ce test repose sur la bi-normalit´e des donn´ees. Attention : la bi-normalit´e est une condition plus forte que la normalit´e
de la premi`ere variable combin´ee a` la normalit´e de la deuxi`eme variable. Pour tester la bi-normalit´e, nous allons faire
appel a` la fonction mshapiro.test du package mvnormtest qui effectue le test de Shapiro-Wilk multidimensionnel.
Comme pour le test unidimensionnel, l’hypoth`ese nulle de ce test est la multinormalit´e des donn´ees tandis que l’hypoth`ese
alternative est la non-multinormalit´e. Avec R, nous pouvons donc effectuer ce test :
>
>
>
>
#Si le package n'est pas encore install´
e
#install.packages("mvnormtest")
library(mvnormtest)
mshapiro.test(rbind(donnees$waistcirc,donnees$DEXfat))
Shapiro-Wilk normality test
data: Z
W = 0.9527, p-value = 0.009369
5
15
●
●
●
0.3
●
●
●
●
10
●
●
●
exemple^2
●
●
●
●
●
● ●
●
0.2
●
●●
●● ●
●
y
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
5
0
−2
0
2
●
●
●
●
● ●●
● ●●
●
0.1
0.0
4
●●
● ●
●●
●
●
● ●
●●
●
●●
●●
●●
● ●
● ●●
●
●
●
●
●
● ● ●
●
●
●
●
● ●
●
●● ● ●
●●
●
●
●●
●
● ●
●● ●
●
●●● ●
●●●
●●
● ●
●
0.0
exemple
0.1
0.2
0.3
x
F IGURE 3 – Une corr´elation nulle n’implique pas qu’il
n’y ait aucun lien entre les variables.
F IGURE 4 – Les donn´ees ne suivent pas une loi binormale
mais l’hypoth`ese de lin´earit´e est v´erifi´ee !
La p-valeur e´ tant inf´erieure au seuil α = 0.05, nous rejetons l’hypoth`ese nulle et acceptons l’hypoth`ese alternative de
non-binormalit´e. En cons´equence, nous ne pouvons pas faire le test sur le cœfficient de corr´elation. Si la p-valeur avait e´ t´e
sup´erieure au seuil α = 0.05 nous aurions pu analyser la p-valeur retourn´ee par la fonction cor.test de R.
Test de lin´earit´e bas´e sur l’anova
Le fait que les donn´ees ne suivent pas une loi bi-normale ne pr´esume en rien de la lin´earit´e des donn´ees, comme le
montre l’exemple suivant, dans lequel nous simulons un mod`ele lin´eaire dont la variable explicative est g´en´er´e grˆace a` une
distribution exponentielle :
>
>
>
>
>
>
set.seed(10)
x<-rexp(100,10)
erreur<-rnorm(100,0,10ˆ-2)
y<-x+erreur
qplot(x,y)
mshapiro.test(rbind(exemple,exemple+erreur))
Shapiro-Wilk normality test
data: Z
W = 0.9985, p-value = 0.5718
Le test de binormalit´e nous ayant conduit a` une impasse, et les donn´ees e´ tant visiblement align´ees, nous allons employer une autre strat´egie. Pour la mettre en œuvre, nous devons d’abord discr´etiser la variable r´eponse waistcirc. Nous
choisissons de faire 5 cat´egories comme suit :
>
>
>
>
>
library(arules)
waist_fact<- discretize(donnees$waistcirc,method="frequency",5)
waist_fact<-as.numeric(waist_fact)
nb_levels<-length(unique(waist_fact))
for(i in 1:nb_levels){
6
indice<-which(i==waist_fact)
waist_fact[indice]<-round(mean(donnees$waistcirc[indice]))
}
> qplot(as.factor(waist_fact),donnees$DEXfat)+
geom_boxplot(outlier.colour = "green",notch=TRUE,
fill = "purple", colour = "#3366FF")
●
●
60
50
donnees$DEXfat
●
●
40
30
●
20
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
10
70
79
86
98
110
as.factor(waist_fact)
F IGURE 5 – Discr´etisation de la variable explicative.
Nous avons cr´ee´ une nouvelle variable, waist fact qui comporte 5 classes. Chaque classe est repr´esent´ee par la valeur moyenne des individus qui la compose. En consid´erant la nouvelle variable tantˆot comme un facteur, tantˆot comme un
vecteur quantitatif, nous allons pouvoir tester la lin´earit´e.
Dans ce cas (c’est-`a-dire quand nous avons plusieurs observations pour chaque niveau de la variable r´eponse), la
r´egression lin´eaire peut eˆ tre vue comme une simplification de l’analyse de la variance. En effet, dans l’analyse de la variance, chaque niveau poss`ede son propre param`etre. Dans la r´egression, les niveaux sont param´etr´es par une droite ; ils ont
donc moins de libert´e. Proc´edons maintenant l’analyse de variance et l’analyse de la variance de la r´egression lin´eaire :
> anov<-aov(donnees$DEXfat˜as.factor(waist_fact))
> anova(anov)
Analysis of Variance Table
Response: donnees$DEXfat
Df Sum Sq Mean Sq F value
Pr(>F)
as.factor(waist_fact) 4
6616 1654.00 56.857 < 2.2e-16 ***
Residuals
66
1920
29.09
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> droite<-lm(donnees$DEXfat˜waist_fact)
> anova(droite)
Analysis of Variance Table
7
Response: donnees$DEXfat
Df Sum Sq Mean Sq F value
Pr(>F)
waist_fact 1 6613.8 6613.8 237.41 < 2.2e-16 ***
Residuals 69 1922.2
27.9
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Les deux tableaux obtenus se ressemblent beaucoup, mais refl`etent deux analyses bien diff´erentes.
La premi`ere analyse est l’analyse de la variance classique, o`u la variable waist fact est consid´er´e comme un facteur
comprenant 5 niveaux. Le nombre de degr´e de libert´e associ´e au facteur est 5 − 1 = 4.
La deuxi`eme analyse o`u waist fact est consid´er´e comme quantitatif. C’est l’analyse de la variance de la r´egression.
Le nombre de degr´e de libert´e correspondant a` la variable explicative est 1.
Les deux tableaux montre que l’anova a un meilleur pouvoir explicatif que la r´egression. La somme des carr´es li´ee au
facteur est plus grande dans le cas de l’anova que de la r´egression. La Figure 6 vous aidera a` comprendre ce que nous venons
d’´ecrire :
> qplot(waist_fact,donnees$DEXfat)+
geom_abline(intercept=coef(droite)[1],slope=coef(droite)[2])+
geom_point(aes(x=c(70,79,86,98,110),
y=unlist(lapply(split(donnees$DEXfat,as.factor(waist_fact)),mean)),
color="red"),lwd=6,alpha=0.8)+
geom_point(aes(x=c(70,79,86,98,110),
y=sort(unique(round(fitted(droite),2))),color="green"),lwd=6,alpha=0.8)+
scale_colour_discrete(name ="",
labels=c("Regression", "ANOVA"))
●
●
60
50
●
●
●
donnees$DEXfat
●
●
40
●
●
●
●
●
●
30
●
20
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
Regression
ANOVA
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
10
70
80
90
100
110
waist_fact
F IGURE 6 – Anova contre r´egression.
´
Ecrivons
maintenant l’´egalit´e de la somme des carr´es pour l’analyse de la variance. Nous avons :
SCT = SCf + SCR ,
8
o`u SCT est la somme des carr´es totale, SCf est la somme des carr´es du facteur (ou encore, somme des carr´es intergroupes) et SCR est la somme des carr´es r´esiduelles (ou encore, somme des carr´es intra-groupes). La quantit´e SCf contient
en r´ealit´e deux e´ l´ements distincts : elle contient la variabilit´e li´ee a` la r´egression lin´eaire, et la variabilit´e qui n’est pas li´ee
a` la r´egression lin´eaire SCl et que nous appelons variabilit´e non-lin´eaire SCnl . Nous pouvons alors r´ee´ crire l’´egalit´e de la
somme des carr´es :
SCT = SCl + SCnl + SCR .
Nous aboutissons au tableau de l’analyse de la variance suivant :
lin´eaire
non-lin´eaire
r´esiduelle
totale
SC
6614
2
1920
8536
ddl
1
3
66
70
CM
6614
0.67
29
Fobs
228
0.03
p-valeur
0
0.99
Les p-valeurs ont e´ t´e obtenues par la commande :
> 1-pf(228,1,66)
[1] 0
> 1-pf(0.03,3,66)
[1] 0.9929393
Pour pouvoir les interpr´eter, il faut v´erifier les conditions fondamentales de l’ANOVA :
> bartlett.test(residuals(anov),as.factor(waist_fact))
Bartlett test of homogeneity of variances
data: residuals(anov) and as.factor(waist_fact)
Bartlett's K-squared = 8.8038, df = 4, p-value = 0.06619
> shapiro.test(residuals(anov))
Shapiro-Wilk normality test
data: residuals(anov)
W = 0.9818, p-value = 0.3921
Les deux p-valeurs e´ tant plus grandes que 0.05, nous d´ecidons de garder les hypoth`eses nulles des deux tests qui sont
la normalit´e des r´esidus et l’homosc´edasticit´e. Nous pouvons maintenant interpr´eter le tableau de l’anova que nous venons
d’´etablir. Les deux p-valeurs correspondent aux deux tests suivants (donn´es dans l’ordre du tableau) :
H0 : il n’existe pas de relation lin´eaire,
contre
H1 : il existe une relation lin´eaire.
et le deuxi`eme test :
9
H0 : le mod`ele s’´ecrit β0 + β1 x1 ou encore il n’existe pas de relation non-lin´eaire,
contre
H1 : le mod`ele ne s’´ecrit pas β0 + β1 x1 ou encore il existe une relation non-lin´eaire.
La p-valeur du premier test est 0, et nous rejetons donc l’hypoth`ese nulle et d´ecidons qu’il existe une partie lin´eaire.
Par ailleurs, la p-valeur du deuxi`eme test est de 0.99, nous gardons donc l’hypoth`ese nulle, et d´ecidons que le mod`ele
s’´ecrit uniquement sous la forme d’une droite affine. Nous pouvons donc, au vu de ces deux tests, d´ecider de la lin´earit´e des
donn´ees.
Le test de Rainbow
Le test de Rainbow permet e´ galement de tester la lin´earit´e des donn´ees. L’id´ee de ce test est de consid´erer que, quand
bien mˆeme les donn´ees ne seraient pas lin´eaires une sous-partie de ces donn´ees peut eˆ tre lin´eaire. Ce test a pour hypoth`ese
nulle la lin´earit´e des donn´ees, contre l’hypoth`ese alternative de non-lin´earit´e des donn´ees. Ce test est fourni par le package
lmtest :
> library(lmtest)
> raintest(lm(DEXfat˜waistcirc,data=donnees))
Rainbow test
data: lm(DEXfat ˜ waistcirc, data = donnees)
Rain = 0.6452, df1 = 36, df2 = 33, p-value = 0.8998
La p-valeur e´ tant sup´erieure a` 0.05, nous d´ecidons de ne pas rejeter l’hypoth`ese nulle, et donc, nous d´ecidons de garder
l’hypoth`ese de lin´earit´e de donn´ees.
2.1.2
Les variables sont mesur´ees sans erreur
Il n’existe pas de mani`ere simple de tester si les variables sont mesur´ees sans erreur. Cependant, nous verrons dans la
partie concernant les diagnostics de la r´egression, que nous v´erifirons si certains points ne jouent pas un trop grand rˆole
dans la r´egression. Ce qu’il faut savoir, c’est que la pr´esence d’erreur de mesure conduit a` une sous-estimation du lien de
lin´earit´e. Nous allons mesurer ce lien de lin´earit´e grˆace a` la corr´elation lin´eaire.
>
>
>
>
>
>
>
x_sans_bruit<-rnorm(100,0,2)
x_bruite_1<-x_sans_bruit+rnorm(100,0,0.5)
x_bruite_2<-x_sans_bruit+rnorm(100,0,1)
x_bruite_3<-x_sans_bruit+rnorm(100,0,1.5)
erreur<-rnorm(100,0,0.1)
y<-x_sans_bruit+erreur
cor(y, x_sans_bruit)
[1] 0.9987345
> cor(y, x_bruite_1)
[1] 0.9681921
> cor(y, x_bruite_2)
[1] 0.88846
> cor(y, x_bruite_3)
10
[1] 0.7543304
La corr´elation diminue lorsque le niveau du bruit (ici, repr´esent´e par l’´ecart type de la loi normale qui augmente)
augmente.
2.1.3
Tester les hypoth`eses implicites
Les hypoth`eses implicites, comme nous l’avons vu, portent sur les termes d’erreur. Pour ce faire, nous devons d’abord
r´ecup´erer dans R le vecteur des r´esidus.
> droite<-lm(DEXfat˜waistcirc,data=donnees)
> residus<-residuals(droite)
En R, lm est la fonction qui permet d’ajuster un mod`ele lin´eaire. La fonction residuals permet de r´ecup´erer les
valeurs r´esiduelles, c’est-`a-dire :
resi = yi − yˆi ,
o`u yˆi = βˆ0 + βˆ1 x1i est la valeur ajust´ee de yi .
Ind´ependance
L’ind´ependance des r´esidus doit eˆ tre, a` ce stade, diagnostiqu´ee en fonction des connaissances que nous avons sur le jeu
de donn´ees. Dans notre cas, nous la supposerons v´erifi´ee. L’analyse des r´esidus dans la partie des diagnostics pourra nous
fournir quelques indications suppl´ementaires.
Cependant, si une forme de d´ependance est soupc¸onn´ee dans les donn´ees, certains tests sp´ecifiques peuvent eˆ tre mis
en place. En particulier, si l’on soupc¸onne une autocorr´elation des erreurs (ce qui peut arriver si les donn´ees sont mesur´ees
temporellement), il existe le test de Breusch-Godfrey, qui peut eˆ tre r´ealis´e a` partir de la commande bgtest du package
lmtest.
Normalit´e
Nous effectuons le test de Shapiro-Wilk, dont l’hypoth`ese nulle est la normalit´e des donn´ees et l’hypoth`ese alternative
la non-normalit´e des donn´ees :
> shapiro.test(residus)
Shapiro-Wilk normality test
data: residus
W = 0.9712, p-value = 0.1016
La p-valeur de ce test e´ tant plus grande que 0.05, nous d´ecidons de ne pas rejeter l’hypoth`ese nulle, et de fait, de d´ecider
que les r´esidus suivent bien une loi normale. Ce choix est fait avec un risque de seconde esp`ece que nous ne calculerons pas
ici.
Homog´en´eit´e des variances
L’homog´en´eit´e des variances peut eˆ tre test´ee grˆace au test de Breusch-Pagan. Il a pour hypoth`ese nulle l’homosc´edasticit´e
et pour hypoth`ese alternative l’h´et´erosc´edasticit´e. Nous utilsons la fonction bptest du package lmtest.
> library(lmtest)
> bptest(droite)
11
studentized Breusch-Pagan test for homoscedasticity
data: droite
BP = 3.1948, df = 1, p-value = 0.07387
La p-valeur e´ tant sup´erieure a` 0.05, nous d´ecidons de garder l’hypoth`ese nulle d’homosc´edasticit´e. Un autre test
e´ quivalent est le test de White, utilisable par la fonction white.test du package bstats ; les hypoth`eses nulle et
alternative sont les mˆeme que pour le test de Breusch-Pagan :
> library(bstats)
> white.test(droite)
White test for constant variance
data:
White = 3.721, df = 2, p-value = 0.1556
La p-valeur e´ tant sup´erieure a` 0.05, nous parvenons a` la mˆeme conclusion et gardons l’hypoth`ese d’homosc´edasticit´e.
Cependant la partie des diagnostics nous permettra e´ galement de v´erifier cette hypoth`ese.
2.2
R´egression lin´eaire et tests
Comme nous l’avons d´ej`a vu, pour appliquer la r´egression lin´eaire nous devons faire :
> droite<-lm(DEXfat˜waistcirc,data=donnees)
La fonction summary permet d’avoir toutes les informations essentielles sur la r´egression :
> summary(droite)
Call:
lm(formula = DEXfat ˜ waistcirc, data = donnees)
Residuals:
Min
1Q
-12.5677 -3.6386
Median
-0.3711
3Q
2.8719
Max
17.7140
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -31.03158
3.67831 -8.436 3.18e-12 ***
waistcirc
0.70740
0.04157 17.017 < 2e-16 ***
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.879 on 69 degrees of freedom
Multiple R-squared: 0.8076,
Adjusted R-squared:
F-statistic: 289.6 on 1 and 69 DF, p-value: < 2.2e-16
0.8048
ˆ 0
Analysons la sortie. Dans la partie cœfficient, nous trouvons deux lignes : l’intercept, qui correspond a` notre beta
ˆ
et waistcirc qui correspond au cœfficient de la variable explicative β1 . La premi`ere colonne, estimate, donne les
estimations ponctuelles pour ces deux estimateurs. Nous pouvons les retrouver grˆace a` la fonction coef :
> coef(droite)
(Intercept)
-31.0315836
waistcirc
0.7073954
12
La colonne Std. Error donne l’´ecart-type pour chacun des estimateurs. Les colonnes suivantes correspondent alors
aux deux tests β0 = 0 contre β0 6= 0 et β1 = 0 contre β1 6= 0. La colonne t value donne la statistique de test, tandis que
la colonne Pr(>|t|) donne la p-valeur associ´ee. Nous avons ensuite, dans l’ordre :
— Une estimation de σ, obtenue par la formule suivante :
> sqrt(sum(residusˆ2)/69)
[1] 4.878985
— Le R2 (et sa version ajust´ee) qui repr´esente la propotion de variabilit´e expliqu´ee par le mod`ele,
— La statistique F de l’analyse de la variance de la r´egression, elle permet de tester la significativit´e globale du mod`ele.
Elle peut e´ galement eˆ tre obtenue de la mani`ere suivante :
> anova(droite)
Analysis of Variance Table
Response: DEXfat
Df Sum Sq Mean Sq F value
Pr(>F)
waistcirc 1 6893.5 6893.5 289.59 < 2.2e-16 ***
Residuals 69 1642.5
23.8
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Comme nous avons v´erifi´e les conditions d’application de la r´egression lin´eaire, nous pouvons faire le test suivant :
H 0 : β1 = 0
contre
H1 : β1 6= 0.
La p-valeur pour ce test e´ tant inf´erieure a` 0.05, nous rejetons l’hypoth`ese nulle et acceptons l’hypoth`ese alternative
β1 6= 0. La variable waistcirc a donc un effet lin´eaire significatif sur la variable DEXfat.
Nous pouvons obtenir des intervalles de confiance pour les param`etres :
> confint(droite)
2.5 %
97.5 %
(Intercept) -38.3696199 -23.6935473
waistcirc
0.6244669
0.7903239
De la mˆeme mani`ere, nous pouvons obtenir des intervalles de confiance pour les valeurs ajust´ees yˆi :
> pred<-predict(droite,interval="confidence")
> head(pred)
ind
ind
ind
ind
ind
ind
1
2
3
4
5
6
fit
39.70795
39.35426
36.87837
19.90088
32.28030
28.03593
lwr
38.14941
37.82324
35.52008
18.17992
31.11191
26.83675
upr
41.26650
40.88528
38.23667
21.62185
33.44870
29.23511
Maintenant, nous pouvons mettre toutes ces informations sur un graphique (Figure 7) obtenu par la commande suivante :
> qplot(donnees$waistcirc,donnees$DEXfat)+
geom_abline(intercept=coef(droite)[1],slope=coef(droite)[2],color="red")+
geom_point(aes(x=donnees$waistcirc,y=fitted(droite)),col="red",lwd=4)+
geom_ribbon(aes(ymin=pred[,2],ymax=pred[,3]),alpha=0.3,fill="green")
13
●
●
60
●
50
●●
donnees$DEXfat
●
●
●●
40
●
●
30
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●
●
● ●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
20
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●● ●
●
●
●
●
●●
●
●
●
●● ●●
●
●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
10
70
80
90
100
110
donnees$waistcirc
F IGURE 7 – La droite de r´egression, avec son intervalle de confiance, et les valeurs ajust´ees.
2.3
Pr´ediction en r´egression lin´eaire
Supposons maintenant que nous disposions de nouveaux individus pour lesquels nous aurions m´esur´e leur tour de taille.
Cr´eeons articiellement un tel vecteur :
> waist_nouv<-c(72.1,90.4,100,111.3)
Nous aimerions savoir, compte tenu de l’analyse de r´egression lin´eaire faite pr´ec´edemment, quelles sont les valeurs
probables que peuvent prendre les mesures de graisse corporelle pour ces individus. Ceci peut eˆ tre fait grˆace a` la fonction
predict avec l’argument newdata. Par ailleurs, l’argument interval permet d’obtenir l’interval de pr´ediction :
> predict(droite,newdata=data.frame(waistcirc=waist_nouv),interval="prediction")
1
2
3
4
2.4
fit
19.97162
32.91696
39.70795
47.70152
lwr
10.08841
23.11215
29.85065
37.70125
upr
29.85484
42.72177
49.56526
57.70180
Diagnostics de la r´egression lin´eaire
Vous avez sans doute not´e que la plupart des tests servant a` v´erifier les conditions d’application du mod`ele ont conduit a`
la non-rejection de l’hypoth`ese nulle. Les d´ecisions que nous avons alors prises ont alors e´ t´e faite avec un risque de seconde
esp`ece que nous ne connaissons pas. Par cons´equent, il n’est pas inutile d’effectuer quelques v´erifications suppl´ementaires.
De plus les diagnostics vont nous permettre de d´etecter des observations aberrantes, et attirer notre attention sur de
possibles probl`emes dans le jeu de donn´ees. C’est donc une e´ tape indispensable.
14
2.4.1
Graphe des r´esidus studentis´es
Nous allons dans ce paragraphe e´ tudier le graphique des r´esidus studentis´es. Les r´esidus, mˆeme sous l’hypoth`ese d’homosc´edasticit´e des erreurs, n’ont pas la mˆeme variance ; ils ne sont donc pas directement comparables. Les r´esidus studentis´es permettent de r´epondre a` cette probl´ematique. Ils sont obtenus de la mani`ere suivante :
> library(MASS)
> residus_stud<-studres(droite)
Nous trac¸ons alors le graphique des r´esidus studentis´ees en fonction des valeurs ajust´ees (Figure 8) :
> qplot(fitted(droite),residus_stud)+
geom_abline(intercept=2,slope=0,color="red",lwd=1.5)+
geom_abline(intercept=-2,slope=0,color="red",lwd=1.5)+
geom_text(aes(label=ifelse((abs(residus_stud)>2),
paste(names(residus_stud), "\n", round(residus_stud,1)),"")), hjust=1.1)
ind 41
●
4.1
4
ind 48
●
2.3
2
residus_stud
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
● ● ●●
●
●
●●
●
0
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
−2
ind 55
●
−2.7
20
30
40
50
fitted(droite)
F IGURE 8 – R´esidus studentis´es contre valeurs ajust´ees.
D’abord, nous pouvons consid´erer que les observations dont les r´esidus studentis´es sont sup´erieures a` 2 en valeur absolue sont suspectes. Dans notre cas nous en avons trois. Nous y reviendrons plus tard.
Une chose importante est d’´etudier la forme du nuage de points, afin de d´etecter des probl`emes de non-lin´earit´e ou
d’h´et´erog´en´eit´e des variances. Pour comparaison, nous produisons les cas les plus classiques dans la Figure 9, et dont le
code est ci-dessous :
>
>
>
>
>
>
library(gridExtra)
gr1<-qplot(1:50,rnorm(50))
gr2<-qplot(1:50,(1:50)*rnorm(50))
gr3<-qplot(1:50,sqrt((1:50))*rnorm(50))
gr4<-qplot(1:50,cos((1:50)*pi/25)+rnorm(50))
grid.arrange(gr1,gr2,gr3,gr4,ncol=2)
15
60
●
●
●
2
●
●
rnorm(50)
●
●● ●
1
●
●
● ●
●
●●
●
●
●
●
0
●
● ●
●
●
●
●
● ●●
●
●
●●
●
●
● ●
●
●
●
● ●
●
●●
−1
(1:50) * rnorm(50)
●
●
●
●
●
● ●
30
●
●
●
●●
●
●● ●● ●
●
●
● ●
●●
0
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●●
●
●
−30
●
●
●
●
●
●
●
●
●
−60
●
−2
0
10
20
30
40
50
●
0
10
20
sqrt((1:50)) * rnorm(50)
10
●
●
●●
5
●
0
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●●
●●●
−5
●
●
●
●
●
●
−10
●
0
10
20
30
40
50
1:50
cos((1:50) * pi/25) + rnorm(50)
1:50
30
40
50
●
3
●
●
2
●●
●
● ●
●
1
●
●
●
●
●
●
●
●
●
●
●
●
● ● ●●
−1
●
●
●
●
0
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
−2
●
●
●
0
10
20
1:50
30
40
50
1:50
F IGURE 9 – Diff´erents cas pour le graphique des r´esidus.
En haut a` droite, c’est la situation ”normale”. Les deux graphiques suivants pr´esentent un cas d’h´et´erosc´edasticit´e, plus
marqu´e en haut a` droite qu’en bas a` gauche. Cele se caract´erise par des formes d’entonnoir ou de diabolo, par exemple. Le
dernier graphique pr´esente un cas de non lin´earit´e.
Dans notre cas, nous pouvons suspecter un l´eger probl`eme d’h´et´etosc´edascit´e.
2.5
Valeurs potentiellement influentes
Pour cela, nous allons calculer ce que nous appelons les hatvalues. Ces valeurs sont li´ees a` la distance qui s´epare
l’observation de la variable explicative avec la moyenne des observation. Il faut les comparer avec la valeur 4/N , avec N le
nombre d’observations. Dans notre cas la valeur de r´ef´erence est : 4/71.
> hv<-hatvalues(droite)
> indice<-1:71
> qplot(indice,hv)+
geom_abline(intercept=0.057,slope=0,col="red",lwd=1.4)+
geom_text(aes(label=ifelse((hv>0.057 &indice!=48),
paste(names(hv), "\n", round(hv,2)),"")), hjust=+1.1)+
geom_text(aes(label=ifelse((hv>0.057 &indice==48),
paste(names(hv), "\n", round(hv,2)),"")), hjust=-0.1)
16
0.08
ind 18
●
0.08
ind 45
●
0.08
ind 47 ind 48
●●
0.07 0.07
ind 46
●
0.06
●
0.06
hv
●
●● ●
●●
●
0.04
●
●
●
●
●
●
●
● ●
●
●
● ●
●
●
●
●
●
●
0.02
●
●
●
●
●
●●
●
●
●
●
●
0
●
●
●
●
●
●
●
●
●
●
●● ●●
●
●
●
●
20
●
●
●
●
●
● ●
●●
40
60
indice
F IGURE 10 – D´etection des valeurs aberrantes.
La Figure 10 indique la pr´esence de 5 individus qui d´epassent le seuil et qui sont donc potentiellement influentes de part
leur e´ loignement de la moyenne des variables explicatives.
2.6
Distance de Cook
La distance de Cook permet de prendre en compte l’´eloignement par rapport a` la moyenne des variables explicatives et
l’´eloignement par rapport a` la droite de r´egression. Les observations dont la distance de Cook est sup´erieure a` 4/N doivent
retenir notre attention.
> cd<-cooks.distance(droite)
> qplot(indice,cd)+
geom_abline(intercept=0.057,slope=0,col="red",lwd=1.4)+
geom_text(aes(label=ifelse((cd>0.057 ),
paste(names(cd), "\n", round(cd,2)),"")), hjust=+1.1)
2.7
Diagnostic final
Nous construisons un graphique permettant de visualiser tous les points r´ep´er´es lors de l’´etape de diagnostic.
> type<-matrix("normal",71,2)
> type[c(14,41,45,48,55),1]<- "Cook or stud.res"
> type[,2]<-rep("",71)
> type[c(18,45,46,47,48),2]<- "hi"
> qplot(data=donnees,waistcirc,DEXfat,geom="line")+
geom_point(aes(color=paste(type[,1],type[,2],sep=" ")),lwd=6)+
geom_smooth(method="lm")+
scale_colour_discrete(name ="",
labels=c("Cook ou stud.res", "Cook ou stud.res et hi","normal","hi"))
17
ind 41
●
0.29
0.3
0.2
cd
ind 48
●
0.18
0.1
ind 45
●
0.06● ●
ind 14
●
0.06
●
0.0
● ●
●● ●●●
●
●
●
●
●●● ●
0
●
● ●
●
●
●
●
●● ● ●
●● ● ● ● ● ●
● ●● ●
● ●
●
20
●
●● ●
●●
●
●
●
●●●●
●
● ● ●●●
40
●
●
●
●
●
60
indice
F IGURE 11 – Distance de Cook.
●
60
●
DEXfat
50
40
30
20
10
●●
●
●
●
●
●
●
●● ●
●
● ●●
●
●●
●
●●
●
●
● ●
●
●● ● ●
●●
●● ●
●
●
●●
●
●
●
●● ●
●
● ●
●● ● ●
●
●
●
●●
● ● ● ●
●
●
●
●
●
70
80
90
100
● Cook ou stud.res
● Cook ou stud.res et hi
● normal
● hi
110
waistcirc
F IGURE 12 – Diagnostic final.
18
© Copyright 2024 ExpyDoc
