Chapitre 2 : Traitement et organisation des données
3.5. Modèles de régressions linéaires : les moindres carrés
3.5.1. But des modèles de régressions linéaires
Evaluer l’influence d’une ou de plusieurs variables indépendantes (quantitatives ou qualitatives) sur une variable dépendante (quantitative continue ou discrète).
Les variables « indépendantes » sont aussi appelées « variables explicatives » La variable « dépendante » est aussi appelée « variable expliquée »
Exemple: Mesurer l’influence des caractéristiques biologiques et les antécédents cliniques des patients sur le niveau de la pression diastolique.
Pour élaborer un tel modèle, distinguons deux cas: l’un où toutes les variables explicatives sont des variables quantitatives (ex: âge, poids, taille, etc..) et l’autre où les variables explicatives incluent à la fois des variables quantitatives et des variables qualitatives (comme le sexe, pratique du sport, fumeur, etc…).
3.5.2. Estimation du modèle dans le cas où toutes les variables explicatives sont quantitatives
myreg1 <- lm(diastolic ~ age + poids +glucose +taille , data
= mydata)
summary(myreg1) # Table de régression
confint(myreg1) # Intervalle de confiance des paramètres
Pour avoir la liste de tous attributs de cette régression, on fait:
attributes(myreg1)
3.5.3. Estimation du modèle dans le cas où les variables
explicatives contiennent des variables qualitatives
113 Lorsque les variables explicatives contiennent des variables catégorielles (ex: sexe, classe_age, imc_cat, tabac, sport), celles-ci doivent être d’abord transformées en variables factors avant d’être intégrées au modèle.
Pour transformer les variables qualitatives en factors on utilise la fonction as.factor() comme suit: +taille+imc_cat+tabac+sport , data = mydata)
summary(myreg2)
NB: Il faut noter que R choisit par défaut la première modalité de la variable catégorielle comme modalité de référence. On peut fixer la modalité de référence en utilisant la fonction relevel(). Ex:
mydata$classe_age <- relevel(mydata$classe_age, ref = "35-45 ans") # les 35-45 ans pris comme reférence.
On peut alors estimer de nouveau le modèle après cette redéfinition.
myreg2 <- lm(diastolic ~ sexe+ classe_age+ poids +glucose +taille +imc_cat + tabac+sport , data = mydata)
summary(myreg2)
3.5.4. Insérer un terme quadratique dans le modèle: la fonction I()
On soupçonne par exemple une relation non linéaire entre l’âge et la pression diastolique. Pour cela, on intègre au modèle un terme quadratique pour capter cette non linéarité. Le terme quadratique correspond à l’âge au carré. Ainsi, on a:
myreg3 <- lm(diastolic ~ sexe+ age+I(age^2)+ poids +glucose +taille+ imc_cat+tabac+sport , data = mydata)
summary(myreg3)
Il est également possible de prendre en compte les effets d’interaction entre les variables. Par exemple l’effet d’interaction entre le niveau de glucose et le tabac. On a:
myreg3b <- lm(diastolic ~ sexe+ age+I(age^2)+ poids +glucose+
tabac+ I(glucose*tabac)+taille+imc_cat +sport , data = mydata) summary(myreg3b)
3.5.5. Sélection automatique des variables explicatives: la fonction
step()
114 La fonction step() permet de sélectionner les variables explicatives pertinentes (sur la base de leur significativité) pour estimer le modèle. L’utilisation de cette fonction est surtout utile lorsque le nombre de variables explicatives est très élevé.
On distingue deux modes de sélection des variables: la « forward selection » et la
« backward selection ».
3.5.5.1. Forward selection
step(lm(diastolic ~1 , data = mydata),diastolic ~ sexe+ age+
poids +glucose +taille+ mal_respire+tabac
+sport+emploi+depression,direction="forward")
Le modèle optimal est :lm(formula = diastolic ~ glucose + taille + age, data = mydata) selectedmodel1<-lm(formula = diastolic ~ glucose + taille + age, data = mydata)
summary(selectedmodel1) 3.5.5.2. Backward selection
step(lm(diastolic ~ sexe+ age+ poids +glucose +taille+
mal_respir+tabac+sport+emploi+depression , data = mydata),direction="backward")
Le modèle optimal est :lm(formula = diastolic ~ glucose + taille + age, data = mydata) selectedmodel2<-lm(formula = diastolic ~ glucose + taille + age, data = mydata)
summary(selectedmodel2)
3.5.6. Test de contraintes linéaires
Soit le modèle linéaire suivant:
myreg4 <- lm(diastolic ~ sexe+ age+ poids +glucose +taille+
tabac+sport , data = mydata) summary(myreg4)
On va tester quelques hypothèses linéaires (formulées sous forme de propositions.
Pour cela, on va utiliser la fonction linearHypothesis() du package car Mise en œuvre des tests
install.packages ("car") library(car)
P1: L'âge n'a aucune influence sur la pression diastolique.
linearHypothesis(myreg4, "age=0") # hypothèse rejetée P2: L'effet total du poids et de la taille sur la pression diastolique vaut 1.5 linearHypothesis(myreg4, "poids+taille=1.5")
P3: Est-ce l'âge et le poids ont le même effet sur la pression diastolique ?
115 linearHypothesis(myreg4, "age=poids")
P4: Est-ce le tabac et le sport ont des effets strictement contraires sur la pression diastolique ?
linearHypothesis(myreg4, "tabac1= -sport1")
3.5.7. Diagnostic des résidus
3.5.7.1. Rappel des hypothèses de base du modèle linéaire Spécification du modèle
𝑦 = 𝛽0+ 𝛽1𝑥1+ 𝛽2𝑥2+ ⋯ + 𝛽𝑘𝑥𝑘+ 𝜀
Les hypothèses qui conditionnent la validité de l’estimation par les MCO sont:
H1: Les erreurs sont de moyenne nulle H2: La variance des erreurs constante
H3: La covariance des erreurs nulle (absence de corrélation des erreurs) H4: Absence de corrélation entre les erreurs et les variables explicatives H5: Les variables explicatives non aléatoires (mesurées sans erreur)
Lorsque l’une de ces hypothèses n’est pas vérifiée, il faut prendre des dispositions méthodologiques pour apporter des corrections.
3.5.7.2. Diagnostic sur la normalité des résidus Soit le modèle:
myreg5 <- lm(diastolic ~ sexe+ age+ poids +glucose +taille+
tabac+sport , data = mydata) summary(myreg5)
Graphiques post-régression:
graphics.off() x.11()
par(mfrow = c(2,2))
plot(myreg5) # Trace 4 graphiques permettant d’analyser le résultat de la régression
3.5.7.3. Histogramme des résidus
resid<-myreg5$residuals # Récupère les résidus de régression m<-mean(resid, na.rm = TRUE) # Calcule la moyenne des résidus std<-sqrt(var(resid, na.rm = TRUE)) # Calcule l'écart-type des résidus
graphics.off();x.11()
116 hist(resid,breaks=seq(min(resid),max(resid),by=(
range(resid)[2]-range(resid)[1])/50), freq = FALSE) # trace l'histogramme
curve(dnorm(x, mean=m, sd=std), col="darkblue", lwd=2, add=TRUE, yaxt="n") # ajoute le courbe de la loi normale.
3.5.7.4. Test de normalité
ks.test(resid, y="pnorm",alternative = "two.sided", exact = NULL)
#Test de shapiro-wilk
shapiro.test(resid)# Shapiro-Wilk Normality Test (stats) 3.5.7.5. Test d’autocorrélation de Durbin-Watson
install.packages("lmtest") library(lmtest)
dwtest(myreg5)
3.5.7.6. Non constance de la variance des erreurs (hétérocédasticité) Test d’hétéroscédasticité de Breush Pagan: le test classique
library(lmtest) bptest(myreg5)
Test d’hétéroscédasticité de Breush Pagan: Non-constant Variance Score Test library(car)
ncvTest(myreg5)
3.5.7.7. Correction de la corrélation des erreurs et de l’hétérocédasticité par les moindres carrés généralisés:
3.5.7.8. Test de multicolinéarité entre les variables explicatives : les Variance Inflation Factors(VIF)
library(car)
vif(myreg5) # affiche les VIFs
117
3.5.8. Prévision à partir du modèle estimé
On peut utiliser les coefficients estimés pour faire une prédiction de la variable dépendante sur un nouvel échantillon. Par exemple, pour prédire la pression diastolique de nouveaux patients dont on connait déjà les caractéristiques biologiques et antécédents cliniques mais pour lesquels la pression diastolique n’a pas encore relevée.
Pour cela, on se sert la fonction predict() et du modèle final estimé.
Exemple:
Soit le modèle final suivant:
finalreg <- lm(diastolic ~ sexe+ age+ poids +glucose +taille+
tabac+sport , data = mydata) summary(finalreg)
Soit un échantillon fictif constitué des dix dernières observations de l’échantillon patientdata:
newsample <- mydata[1:10,]
Prévision sur le nouvel échantillon
predval<-predict(finalreg, newdata = newsample) # Renvoie les valeurs prédites pour les 10 observations
predval # vecteur des pressions diastoliques prévues