Régression linéaire simple : quand les hypothèses ne sont pas satisfaites
Pour faire suite au précédent post sur la régression linéaire simple dans cet article, je vais vous parler des alternatives possibles à cette méthode, lorsqu’une, ou plusieurs de ses hypothèses de validité, ne sont pas satisfaites.
Table des matières
Introduction
Pour rappel, la régression linéaire simple est une méthode statistique très classique qui est employée pour évaluer si deux variables numériques continues sont significativement liées, en faisant l’hypothèse que la forme de leur relation est linéaire.
La régression linéaire simple permet d’estimer les paramètres de la droite liant la variable réponse à la variable prédictive, mais elle permet également d’évaluer si cette relation est significative ou non. Pour cela, un test T est employé pour évaluer si la pente est significativement différente de 0 ou non. Cependant, le résultat de ce test T n’est valide que si les quatre hypothèses suivantes sont satisfaites :
1. La relation entre les deux variables numériques continues est réellement de forme linéaire (au moins grossièrement).
2. Les résidus de la régression linéaire sont indépendants.
3. Les résidus de la régression linéaire sont distribués selon une loi Normale de moyenne 0 (au moins grossièrement).
4. Les résidus de la régression linéaire sont distribués de façon homogène.
Comme décrit dans le premier article (faire un lien), la linéarité peut être évaluée graphiquement grâce à la fonction “scatterplot” du package car. L’hypothèse d’indépendance, quand à elle, peut être jugée par un “lag plot”, ainsi qu’en employant le test de Durbin-Watson. L’hypothèse de normalité peut être appréciée par la réalisation d’un “QQplot”, couplé à un test de Shapiro-Wilk. Enfin, l’hypothèse d’homogénéité des résidus peut être testée par la réalisation d’un “residuals vs fitted plot“, complété par un test de Breush-Pagan.
J’associe très souvent une méthode d’évaluation visuelle avec la réalisation d’un test statistique pour évaluer ces hypothèses, car les deux sont informatifs. Le seuil de significativité, généralement fixé à 0.05, ne doit pas être appliqué strictement, car certains tests, comme celui de Shapiro-Wilk par exemple, sont peu robustes, ou sensibles aux outliers. La p-value du test statistique doit plutôt être envisagée comme un garde fou. De l’autre coté, les évaluations visuelles apportent une impression générale, globale, qui doit également être prise en
considération.
Quand on a beaucoup d’expérience, on peut se baser uniquement sur les représentations visuelles pour évaluer les différentes hypothèses, mais lorsqu’on débute, je sais qu’il est très difficile de différencier un écart acceptable d’un écart trop important. C’est aussi pour acquérir de l’expérience que je vous conseille de faire cette double évaluation (visuelle et par les tests statistiques) des hypothèses.
Et au final, est ce que c’est grave lorsque les hypothèses ne sont pas satisfaites ? Est ce que cela à un impact sur les résultats ? De mon point de vu, il est très difficile de répondre à cette question, car elle dépend de l’hypothèse concernée, de l’amplitude de l’écart, et enfin du nombre de données.
Sans que je puisse apporter de références pour cela, je dirai qu’un léger défaut de linéarité n’est pas très important, que la non indépendance des résidus est rédhibitoire. De même un défaut de normalité est moins grave qu’un défaut d’homogénéité.
Quand l' hypothèse de linéarité n'est pas satisfaite
Comme expliqué précédemment, l’hypothèse de linéarité peut s’évaluer à l’aide de la fonction “scatterplot” du package “car”.
Pour reprendre l’exemple de l’article précédent :
library(car)
scatterplot(prestige~education, data=Prestige) 
Ici l’hypothèse de linéarité est satisfaite
scatterplot(prestige~income, data=Prestige) 
Là, en revanche, elle est plus problématique. Dans le cas d’un défaut de linéarité, il existe deux alternatives : abandonner la linéarité, ou essayer de l’améliorer.
Abandonner la linéarité et utiliser la corrélation de Spearman
La première solution consiste à abandonner la relation linéaire, et à se satisfaire d’une relation monotone entre les deux variables numériques continues. Et par conséquent, d’utiliser un test de corrélation de Spearman pour évaluer s’il existe une dépendance monotone entre les deux variables. Pour cela on utilise la fonction “cor.test” en employant l’option “spearman” dans l’argument “method”.
cor.test(Prestige$prestige,Prestige$income, data=Prestige, method="spearman")
##
## Spearman's rank correlation rho
##
## data: Prestige$prestige and Prestige$income
## S = 43675, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.7530408 L’avantage de la corrélation de Spearman est qu’elle ne requiert que l’indépendance des données.
Ainsi, soit le protocole expérimental est connu et l’indépendance des données peut être établi, soit le protocole n’est pas connu et dans ce cas là, il est préférable de vérifier cette hypothèse. Le plus simple est alors d’ajuster le modèle linéaire simple et d’évaluer l’auto-corrélation des résidus par un “lag plot” avec la fonction “acf” du package “car” (voir plus bas).
Améliorer la linéarité en utilisant une transformation
La deuxième solution consiste à utiliser une transformation de la variable prédictive afin d’améliorer la linéarité. Puis d’utiliser cette variable prédictive transformée lors de la modélisation. Les autres hypothèses doivent ensuite être évaluées.
A mon sens, cette alternative ne doit être envisagée que lorsque l’évaluation du lien linéaire est indispensable.
Les principales transformations utilisées sont le log et la racine carrée. Ces transformations ont également tendance à améliorer la normalité et l’homogénéité des résidus. En pratique la fonction “log10” (logarithme en base 10) ou la fonction “log” (logarithme néperien) peuvent être employées. Néanmoins, je vous conseille d’utiliser la fonction “log1p” qui applique la transformation log(1+x). Cela vous évitera d’aboutir à des valeurs négatives lorsque la valeur de x est proche de 0. (log1p(0.007) = 0.007 alors que log(0.07)=-4.96) Ceci n’est pas forcément embêtant ici, mais peut l’être lorsque la transformation log est appliquée sur une variable réponse.
# evaluation de la linéarité aprés transformation de la variable prédictive
scatterplot(prestige~log1p(income), data=Prestige)
Pour la racine carrée, il suffit d’utiliser la fonction “sqrt”.
scatterplot(prestige~sqrt(income), data=Prestige) 
Ici, la transformation par la racine carrée fonctionne bien.
Dans un second temps, on réalise la régression linéaire avec cette transformation.
mod_sqrt <- lm(prestige~sqrt(income), data=Prestige) summary(mod_sqrt) ## ## Call: ## lm(formula = prestige ~ sqrt(income), data = Prestige) ## ## Residuals: ## Min 1Q Median 3Q Max ## -23.786 -7.897 -2.721 8.075 32.024 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.17653 4.03586 0.539 0.591
## sqrt(income) 0.56387 0.04895 11.519 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.33 on 100 degrees of freedom
## Multiple R-squared: 0.5703, Adjusted R-squared: 0.566
## F-statistic: 132.7 on 1 and 100 DF, p-value: < 2.2e-16 Quand l'hypothèse d'indépendance des résidus n'est pas satisfaite
Comme évoqué plus haut, l’indépendance des résidus peut s’évaluer avec un “lag plot” sur les résidus du modèle de régression linéaire. On peut également utiliser le test de Durbin Watson, mais celui-ci n’évalue que le lag1.
Exemple de non-indépendance des résidus
Un cas classique de non-indépendance des résidus est rencontré lorsque la variable prédictive est une variable temporelle.Dans cette situation, les résidus sont dits “auto-corrélés“.
Pour illustrer cette situation, nous allons utiliser une partie du jeu de données “airquality” du package “dataset” (installé par défaut à chaque session R), en le réduisant aux données du mois numéro 9.
library(tidyverse)
aq9 <- airquality %>%
filter(Month==9)
head(aq9)
## Ozone Solar.R Wind Temp Month Day
## 1 96 167 6.9 91 9 1
## 2 78 197 5.1 92 9 2
## 3 73 183 2.8 93 9 3
## 4 91 189 4.6 93 9 4
## 5 47 95 7.4 87 9 5
## 6 32 92 15.5 84 9 6 Si on ajuste le modèle de régression de la variable “Ozone” en fonction de la variable “Day” (qui va de 1 à 30), on s’aperçoit que les résidus sont auto-corrélés.
# régression linéaire simple
aq9.lm <- lm(Ozone~Day, data=aq9)
# lag plot
acf(residuals(aq9.lm), main="aq9.lm") 
Le lag plot montre clairement que les résidus de lag 1 sont auto-corrélés, puisque le coefficient de corrélation dépasse la borne supérieure de l’intervalle de confiance du coefficient de corrélation de valeur nulle.
durbinWatsonTest (aq9.lm)
## lag Autocorrelation D-W Statistic p-value
## 1 0.4771818 0.8561456 0
## Alternative hypothesis: rho != 0 Le test de Durbin-Watson confirme l’auto-corrélation de lag1, puisque la pvalue est égale à 0.
Cela n’est pas étonnant, car on peut légitimement penser que le niveau d’Ozone du jour j+1 dépend de celui observé au temps précédent, c’est-à-dire au jour j.
Prise en compte de l'auto-corrélation par ajout d'une structure de corrélation
Dans cette situation, il est nécessaire d’ajouter une structure corrélation des résidus, dans le modèle de régression, afin de modéliser cette dépendance. Concrètement, il s’agit d’ajouter une matrice de variance co-variance de type auto-regressive 1 (AR-1). Pour cela, il est nécessaire d’utiliser la fonction “gls” du package “nlme”, et d’ajouter la structure corAR1(form=~Day) dans l’argument “correlation”.
library(nlme)
aq9.gls <- gls(Ozone~Day, data=aq9, correlation=corAR1(form=~Day),na.action=na.omit) Afin de vérifier que l’ajout de cette structure de corrélation améliore la qualité du modèle, on peut comparer les critères d’Akaiké (AIC). Mais comme le modèle avec structure de corrélation a été ajusté avec la fonction gls, il faut refaire tourner un modèle sans structure, avec la fonction “gls” également.
# modèle sans structure de corrélation
aq9.gls0 <- gls(Ozone~Day, data=aq9,na.action=na.omit)
AIC(aq9.gls,aq9.gls0)
## df AIC
## aq9.gls 4 241.1113
## aq9.gls0 3 251.1822 L’AIC du modèle avec structure de variance est nettement plus faible, ce qui montre que sa qualité est meilleure.
Lorsqu’on regarde les sorties de ces deux modèles, on s’aperçoit que le coefficient de la pente a un peu bougé : -1.8 pour le modèle sans structure contre -2.1 pour le modèle avec la structure de corrélation AR1.
#modèle sans structure de corrélation
summary(aq9.gls0)
## Generalized least squares fit by REML
## Model: Ozone ~ Day
## Data: aq9
## AIC BIC logLik
## 251.1822 255.0697 -122.5911
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 59.12181 6.982162 8.467551 0e+00
## Day -1.83227 0.402514 -4.552052 1e-04
##
## Correlation:
## (Intr)
## Day -0.871
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.4219861 -0.8721532 -0.1600917 0.8569013 2.1201810
##
## Residual standard error: 18.49241
## Degrees of freedom: 29 total; 27 residual
#modèle avec structure de corrélation
summary(aq9.gls)
## Generalized least squares fit by REML
## Model: Ozone ~ Day
## Data: aq9
## AIC BIC logLik
## 241.1113 246.2946 -116.5556
##
## Correlation Structure: ARMA(1,0)
## Formula: ~Day
## Parameter estimate(s):
## Phi1
## 0.7366436
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 68.03627 18.170194 3.744389 0.0009
## Day -2.10930 0.978811 -2.154961 0.0402
##
## Correlation:
## (Intr)
## Day -0.834
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.4300307 -0.9337497 -0.4233752 0.4853264 1.3496458
##
## Residual standard error: 23.26605
## Degrees of freedom: 29 total; 27 residual En revanche, l’erreur standard de la pente varie énormément, puisqu’elle est multipliée par un facteur environ 2.5 : elle est égale à 0.4 dans le modèle sans structure de corrélation, alors qu’elle est égale à 0.98 dans le modèle prenant en compte l’auo-corrélation des résidus. L’influence de cette augmentation de l’erreur standard sur la valeur de la statistique T est très importante, puisqu’elle passe de -4.55 dans le modèle sans prise en compte de l’auto-corrélation, à -2.15 lorsque celle-ci est modélisée. Avec évidemment une répercussion sur la p-value qui passe de 1e-04 à 0.04.
De manière générale, lorsque l’auto-corrélation des résidus n’est pas prise en compte, la régression linéaire conduit à des faux positifs, c’est-à-dire à conclure que la pente est significativement différente de 0, alors qu’en réalité elle ne l’est pas.
Autre précision importante, lorsque les résidus sont auto-corrélés, la corrélation de Spearman ne peut pas être employée, car cette méthode nécessite l’indépendance des données (c’est sa seule hypothèse).
Il existe d’autres structures de corrélation, et vous trouverez davantage d’informations sur celles-ci dans le livre “Mixed Effects Models and Extensions in Ecology with R” d’Alain F. Zuur et Elena N. Ieno,
Ainsi que dans le livre “Mixed-Effects Models in S and S-Plus” de Bates et Pinheiro.
Il existe encore d‘autres types de dépendance, comme la dépendance spatiale. L’évaluation se fait par un variogramme. Et dans ce cas, la modélisation de la structure de corrélation est faite en fonction des coordonnées de latitude et longitude des points, en employant des structures de type “corExp”, “corGaus”, “corSpher” ou encore “corLin ” par exemple.
Les structures de corrélation prises en charge par le package “nlme” sont décrites dans la page d’aide de la fonction corClasses.
?corClasses Quand la normalité n'est pas satisfaite
L’hypothèse de normalité des résidus est sans doute l’hypothèse la moins importante pour la régression linéaire simple, car le modèle linéaire est sensé être robuste. Néanmoins, je vous conseille tout de même de l’évaluer visuellement par un QQplot et par un test de Shapiro-Wilk. Si le défaut de normalité est marqué sur le QQ plot, c’est à dire si les résidus sont distribués de façon incurvés, ou si vous avez peu de données (entre 10 et 20) et qu’il y a un ou plusieurs outliers, je vous conseille de vous diriger vers une solution alternative. Il y en a trois, dans l’ordre de mes préférences :
1. Renoncer à la linéarité, se contenter de la monotonie, et utiliser le test du coefficient de corrélation de Spearman, puisque cette méthode ne nécessite que l’hypothèse d’indépendance.
2. Utiliser un test de permutation pour évaluer l’égalité à 0 de la pente puisque les tests de permutations ne nécessitent, eux aussi, que l’hypothèse d’indépendance
3. Appliquer une transformation log, racine carrée ou de type Box Cox sur la réponse afin d’améliorer la normalité des résidus, et refaire tourner le modèle de régression linéaire en appliquant la transformation.
Illustrons ce cas de figure avec le jeu de données suivant :
mydata <-structure(list(x = c(19.79, 8.56, 3.2, 2.33, 5.63, 16.05, 7.46,
19.47, 4.15, 9.72, 4.26, 5.4, 15.68, 2.83, 9.62, 2.61, 11.65,
1.17, 19.73, 7.02), y = c(128.58, 22.57, 12.36, 11.12, 16.89,
37.76, 28.27, 46.22, 16.72, 97.04, 18.16, 175.79, 119.64, 15.08,
70.96, 15.7, 29.01, 9.05, 47.16, 31.01)), class = "data.frame", row.names = c(NA, -20L))
library(ggplot2)
ggplot(mydata, aes(x=x, y=y))+
geom_point()+
theme_classic() 
Lorsqu’on réalise la régression linéaire simple et qu’on évalue l’hypothèse de normalité, on s’aperçoit qu’il existe un défaut manifeste, avec la présence de 4 ou 5 outliers.
mod <- lm(y~x, data=mydata)
plot(mod,2) 
shapiro.test(residuals(mod))
##
## Shapiro-Wilk normality test
##
## data: residuals(mod)
## W = 0.72095, p-value = 7.1e-05 Abandonner la linéarité et utiliser la corrélation de Spearman
cor.test(y,x, data=mydata, method="spearman")
##
## Spearman's rank correlation rho
##
## data: y and x
## S = 258, p-value = 1.985e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.806015 Ici, la p-value est largement inférieure au seuil de significativité généralement employé de 0.05. On conclut donc à la dépendance monotone significative entre les variables x et y.
Utiliser un test de permutation
Ici, la p-value est largement inférieure au seuil de significativité généralement employé de 0.05. On conclut donc à la dépendance monotone significative entre les variables x et y.
library(lmPerm)
set.seed(220718)
modp <- lmp(y~x, data=mydata)
summary(modp)
##
## Call:
## lmp(formula = y ~ x, data = mydata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -36.257 -21.077 -14.520 -2.688 139.567
##
## Coefficients:
## Estimate Iter Pr(Prob)
## x 3.287 2502 0.0388 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 43.34 on 18 degrees of freedom
## Multiple R-Squared: 0.1908, Adjusted R-squared: 0.1459
## F-statistic: 4.245 on 1 and 18 DF, p-value: 0.05412 Si on compare avec les résultats de la régression linéaire simple classique, on s’aperçoit que les p-values sont plutôt proches.
summary(mod)
##
## Call:
## lm(formula = y ~ x, data = mydata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -36.257 -21.077 -14.520 -2.688 139.567
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.471 17.082 1.081 0.2938
## x 3.287 1.596 2.060 0.0541 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 43.34 on 18 degrees of freedom
## Multiple R-squared: 0.1908, Adjusted R-squared: 0.1459
## F-statistic: 4.245 on 1 and 18 DF, p-value: 0.05412 Ce qui va dans le sens de la robustesse de la régression linéaire vis-à-vis d’un défaut de normalité des résidus. Néanmoins, il me semble prudent de le vérifier par un test de permutation, quitte à garder les résultats de la régression s’il ya peu de différence.
Utiliser une transformation sur la variable réponse
Transformation log (fonction log1p)
mod2 <- lm(log1p(y)~x, data=mydata)
plot(mod2,2) 
shapiro.test(residuals(mod2))
##
## Shapiro-Wilk normality test
##
## data: residuals(mod2)
## W = 0.80832, p-value = 0.001153 Ici, la transformation log n’est pas satisfaisante.
Transformation racine carrée (fonction sqrt)
mod3 <- lm(sqrt(y)~x, data=mydata)
plot(mod3,2) 
shapiro.test(residuals(mod3))
##
## Shapiro-Wilk normality test
##
## data: residuals(mod3)
## W = 0.76186, p-value = 0.0002459 La transformation racine carrée n’est pas non plus satisfaisante.
Transformation Box Cox (fonction powerTransform)
La transformation Box-Cox est définie par :
\[
T_{BC}(y,\lambda) = y^{\lambda} \left\{ \begin{array}{rcl}
{\frac{y^{\lambda}-1}{\lambda}} & \mbox{when} & \lambda\neq 0 \\
log(y) & \mbox{when} & \lambda = 0\\
\end{array}\right.
\]
Il s’agit d’une généralisation de la transformation logarithmique. Lorsque le paramètre lambda est égal à 0, alors la transformation Box-Cox revient à une transformation log (népérien).
Le paramètre Lambda est estimé par maximisation de la log vraisemblance. En pratique, on peut réaliser une transformation Box-Cox de la variable réponse, en utilisant la fonction “powerTransform” du package “car”. Le principe est de :
- Estimer le paramètre lambda,
- Créer la variable réponse transformée,
- Refaire tourner le modèle de régression avec la réponse transformé
library(car)
# estimation du paramètre lambda
summary(p1 <-powerTransform(mod))
## bcPower Transformation to Normality
## Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## Y1 -0.5401 -1 -1.026 -0.0541
##
## Likelihood ratio test that transformation parameter is equal to 0
## (log transformation)
## LRT df pval
## LR test, lambda = (0) 4.565952 1 0.032614
##
## Likelihood ratio test that no transformation is needed
## LRT df pval
## LR test, lambda = (1) 31.64622 1 1.8497e-08 Ici, lambda est estimé à -0.5401. Cependant, en pratique, on utilise généralement une valeur arrondie, c’est celle qui est rapportée dans la colonne “Rounded Pwr” ; ici -1. Les deux autres colonnes sont les bornes inférieure et supérieure de l’intervalle de confiance à 95% du paramètre
lambda.
Les lignes en dessous fournissent des informations pour les deux valeurs entières de lambda encadrant le lambda optimal. Lorsque la p-value est inférieure à 0.05, cela signifie que les valeurs de lambda ne sont pas satisfaisantes. Le seuil de 0.05 ne doit pas être appliqué strictement, c’est un repère.
Ici, aucune ne semble satisfaisante. Néanmoins, quand on applique la transformation, l’hypothèse de normalité devient satisfaisante.
# création d'un nouveau jeu de données contenant la variable réponse transformée
mydata2 <- transform(mydata, y2=bcPower(y, coef(p1)),
y2round=bcPower(y, coef(p1, round=TRUE)))
head(mydata2)
## x y y2 y2round
## 1 19.79 128.58 1.717246 0.9922227
## 2 8.56 22.57 1.507653 0.9556934
## 3 3.20 12.36 1.375445 0.9190939
## 4 2.33 11.12 1.347463 0.9100719
## 5 5.63 16.89 1.449348 0.9407934
## 6 16.05 37.76 1.591130 0.9735169
# on refait tourner le modèle avec la variéble réponse transformé avec un lambda arrondi
mod4 <- lm(y2round~x, data=mydata2)
# on visualise la normalité
plot(mod4,2) 
# test de Shapiro-Wilk
shapiro.test(residuals(mod4))
##
## Shapiro-Wilk normality test
##
## data: residuals(mod4)
## W = 0.97177, p-value = 0.7916 On voit ici, que la transformation Box-Cox permet de bien améliorer la normalité des résidus.
Remarque : Lorsque la variable y n’est pas strictement positive, il est nécessaire de rajouter une constante dans la fonction “bcPower” afin que toutes les valeurs soient positives, comme ceci, par exemple en ajoutant 1.
y2=bcPower(y + 1, coef(p1)) De mon point de vu, en cas de rejet de l’hypothèse de normalité, il est préférable d’employer une corrélation de Spearman, ou un test de permutation, car cela permet de continuer à travailler sur les données brutes.
A noter que les transformations log, racine carrée, ou Box-Cox améliore également l’homogénéité des résidus.
Quand l'homogénéité n'est pas satisfaite
L’hypothèse d’homogénéité des résidus est importante, car lorsqu’elle n’est pas satisfaite, l’erreur standard de la pente est également sous estimée. Comme ce paramètre est le dénominateur de la statistique t du test de l’égalité à 0 de la pente, lorsque l’erreur standard de la pente est sous estimée, alors la statistique t est sur-estimée. Et par conséquent la p-value est plus faible qu’elle ne devrait.
$ t\;value= \frac{\widehat{pente}}{\widehat{se_{pente}}} $
L’hypothèse d’homogénéité des résidus peut être évaluée par un “residuals vs fitted” plot, en utilisant la racine carrée des résidus standardisés, ainsi que par le test de Breush-Pagan. Lorsque cette hypothèse est clairement rejetée, une des alternatives suivantes peut être employée :
- Renoncer à la linéarité, se contenter de la monotonie, et utiliser le test du coefficient de corrélation de Spearman, puisque cette méthode ne nécessite que l’hypothèse d’indépendance.
- Utiliser un test de permutation pour évaluer l’égalité à 0 de la pente puisque les tests de permutations ne nécessitent, eux aussi, que l’hypothèse d’indépendance.
- Appliquer une transformation log, racine carrée ou de type Box Cox sur la réponse afin d’améliorer l’homogénéité des résidus, et refaire tourner le modèle de régression linéaire.
- Utiliser des estimateurs sandwich de la matrice de variance-covariance des paramètres du modèle
Les trois premières alternatives ont été présentées dans le chapitre précédent car elles sont aussi des alternatives à la régression linéaire simple en cas de rejet de l’hypothèse de normalité des résidus. Nous allons donc nous focaliser sur l’emploi des estimateurs sandwich.
Utiliser des estimateurs sandwich pour estimer l'erreur standard des paramètres de la régression
Comme expliqué précédemment, en cas d’hétérogénéité des résidus, l’erreur standard des résidus est sous-estimée. Il s’agit donc, d’utiliser une méthode de calcul plus robuste, qui va conduire à une meilleure estimation de l’erreur standard, avec une valeur plus élevée.
En pratique, il s’agit d’estimer la matrice de variance-co-variance des paramètres du modèle de régression avec cette méthode robuste. Puis d’utiliser celle-ci pour calculer la statistique des tests d’égalité à 0 des paramètres, et donc pour dériver les p-values.
Data
ydf <- structure(list(var_rep = c(0.13, 0.15, 0.05, 0.11, 0.1, 0.05,
0.1, 0.02, 0.14, 0.07, 0.07, 0.12, 0.12, 0.33, 0.14, 0.2, 0.16,
0.64, 0.33, 1.11, 0.26, 0.36, 0.39, 0.59, 0.94), var_pred = c(4.36,
3.91, 3.98, 4.24, 3.84, 4.32, 4.65, 4.89, 4.36, 5.1, 4.91, 6.19,
8.05, 4.15, 6.1, 5.32, 5.39, 6.67, 5.42, 10.42, 7.95, 5.99, 5.92,
7.2, 8.77)), class = "data.frame", row.names = c(NA, -25L))
head(mydf)
## var_rep var_pred
## 1 0.13 4.36
## 2 0.15 3.91
## 3 0.05 3.98
## 4 0.11 4.24
## 5 0.10 3.84
## 6 0.05 4.32 Dans l’exemple ci-dessous, l’hypothèse d’homogénéité des résidus est clairement rejetée par la méthode visuelle, puisque celle-ci met en évidence une augmentation de la dispersion des résidus lorsque les fitted values augmentent.
mydf.lm1<-lm(var_rep~var_pred, data=mydf)
plot(mydf.lm1,3) 
Le test de Breush-Pagan rejette également l’hypothèse d’homogénéité des résidus.
ncvTest(mydf.lm1)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 10.43117 Df = 1 p = 0.001239063 Mise en application
Le package “car” contient la fonction “hccm” qui fait appel aux estimateurs sandwich pour calculer la matrice de variance-covariance des paramètres du modèle de régression. On peut facilement comparer les deux types d’estimations avec les commandes suivantes :
#matrice de variance covariance classique biaisée
vcov(mydf.lm1)
## (Intercept) var_pred
## (Intercept) 0.016443387 -0.0026660939
## var_pred -0.002666094 0.0004690524
#matrice de variance covariance de type sandwich, non biaisée
hccm(mydf.lm1)
## (Intercept) var_pred
## (Intercept) 0.034935509 -0.006921778
## var_pred -0.006921778 0.001402946 On peut voir que la variance de la pente passe de 0.00047 à 0.0014, elle a quasiment était multiplié par 3.
Pour obtenir les résultats habituels de la régression, il est nécessaire d’utiliser la fonction “coeftest” du package “lmtest”, comme ceci :
library(lmtest)
coeftest(mydf.lm1, vcov = hccm)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.463146 0.186910 -2.4779 0.020991 *
## var_pred 0.128491 0.037456 3.4305 0.002283 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 A comparer avec la sortie de la fonction summary :
summary(mydf.lm1)
##
## Call:
## lm(formula = var_rep ~ var_pred, data = mydf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.45121 -0.09775 0.02834 0.09672 0.27628
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.46315 0.12823 -3.612 0.00147 **
## var_pred 0.12849 0.02166 5.933 4.77e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1795 on 23 degrees of freedom
## Multiple R-squared: 0.6048, Adjusted R-squared: 0.5876
## F-statistic: 35.2 on 1 and 23 DF, p-value: 4.766e-06 On peut voir que la statistique t de la pente était de 5.9, alors qu’elle n’est plus que de 3.43 lorsque l’hétérogénéité des résidus est prise en compte.
Remarque : En cas d’hétérogénéité des résidus, il est également possible de modéliser leur variance en ajoutant une structure de variance. Cette alternative me semble un peu complexe à mettre en oeuvre alors qu’il existe d’autres solutions plus simples, comme nous venons de le voir.
Néanmoins, pour ceux qui souhaitent se pencher sur la modélisation de la variance, cette méthode est très bien décrite dans le livre “Mixed Effects Models and Extensions in Ecology with R” d’Alain F. Zuur et Elena N. Ieno.
Pour aller plus loin
Pour avoir une vision plus large de cette question de la validité des hypothèses, vous pouvez lire cet article de Lise Vaudor , pour qui le respect de la normalité et de l’homogénéité ne doit pas être pris au pied de la lettre, et cet article de Guillaume Rousselet, qui lui au contraire, se bat pour promouvoir des tests statistiques plus robustes.
Enfin, pour plus d’informations sur les transformations, et plus généralement, sur tout ce qui touche à la régression linéaire, je vous recommande le livre ” An R Companion to Applied Regression” de John Fox et Sanford Weisberg.
Conclusion
Comme je l’ai dit en introduction, de mon point de vu, les hypothèses les plus importantes sont celles de l’indépendance et l’homogénéité des résidus car elles sont une forte influence sur la statistique du test de la pente. L’hypothèse de linéarité est également importante mais je pense qu’une certaine flexibilité peut être conservée, de même pour l’hypothèse de normalité.
J’espère que cet article vous sera utile, et qu’il vous permettra de réaliser vos régressions linéaires simples dans les meilleurs conditions possibles, et que vous saurez quoi faire, et comment, lorsqu’une ou plusieurs hypothèses ne sont pas satisfaites.
Si cet article vous a plu, ou vous a été utile, et si vous le souhaitez, vous pouvez soutenir ce blog en faisant un don sur sa page Tipeee
Note : Je touche une petite commission (entre 3 et 6%) si vous passez par les liens Amazon de cet article pour acheter les livres mentionnés. cela m’aide à entretenir ce blog, merci si vous le faites ! 
Credits photos : geralt
Merci pour cet enseignement de qualité.
Bonjour,
Merci beaucoup pour vos tutos !
Pouvez-vous m’éclairer sur un point, quand on transforme les données, comment pouvons nous interpréter l’équation de la droite ? Cela n’a plus vraiment de sens finalement non ?
Merci 🙂
Bonjour,
les méthodes que vous présenter pour analyser le respect des hypothèses s’appliquent-ils de la même façon aux modèles de régression linéaire multiple?
Merci,
Bonjour Nicolas,
oui les hypothèses s’appliquent aussi aux modèles de régressions linéaires multiples, car elles concernent tous les modèles linéaires.
Bonne continuation.
Bonjour Madame,
Je vous remercie vivement pour vos explications et vos codes qui m’ont éclairé plus d’une fois.
Je me heurte cette fois ci à un problème que j’ai dû mal à identifier, je dispose d’un jeu de donnée avec 33 années et différentes variables quantitatives (une valeur par année pour chacune des variables). Etant donné qu’il y a une auto corrélation des résidus, j’ai ajouté une matrice de variance covariance de type auto régressive 1, comme vous le recommandiez.
Cependant lorsque je fais tourner le code suivant :
reg.s.h.gls <- gls(count~year, data=mhwSub, correlation=corAR1(form=~year),na.action=na.omit)
summary(reg.s.h.gls)
Cela ne me donne aucun coefficient dans la sortie du modèle :
Generalized least squares fit by REML
Model: count ~ year
Data: mhwSub
Correlation Structure: AR(1)
Formula: ~year
Parameter estimate(s):
Phi
0.4497734
Coefficients:
Correlation:
(Intr)
year -1
Standardized residuals:
Min Q1 Med Q3 Max
-1.4590774 -0.8226785 -0.2084742 0.6149282 2.0355680
Avez-vous une idée de ce qui pourrait causer ce problème ? Existe-t-il une autre fonction que gls ?
Je vous remercie,
Cordialement,
Bonjour,
heu….là non comme ça je ne vois pas trop….
Vous trouverez peut être des réponses dans le livre d’Alan Zuur : https://www.amazon.fr/Effects-Models-Extensions-Ecology-Author/dp/B007S7GBTM/ref=sr_1_1?__mk_fr_FR=ÅMÅŽÕÑ&keywords=alain zuur&qid=1578775242&sr=8-1
Bonne continuation
Bonjour et merci beaucoup pour vos articles, toujours très pédagogiques.
Je me pose une question: lorsque l’on fait de la sélection de variables ou de la sélection de modèles, en régression linéaire multiple, est-on est censé tester les hypothèses sur le modèle global incluant toute les variables avant sélection, sur le meilleur modèle obtenu après sélection, ou bien avant et après sélection?
C’est probablement une question stupide mais je n’arrive pas à trouver la réponse!
Bonjour Nicolas,
je suis assez persuadée que c’est loin d’être une question stupide. En pratique on vérifie les hypothèses sur le modèle complet (global), et rien n’empêche de jeter un coup d’oeil (rapide, et en croisant les doigts pour que tout soit OK 😉 ) sur le modèle final.
Bonne continuation.
Bonjour, je découvre votre site qui est super, merci pour tous ces articles !
En cas de relation non linéaire, que pensez vous d’une régression quadratique (y = ax bx² c) au lieu de transformer y ou x ? ça se fait assez facilement dans R, et on a l’avantage de garder les valeurs de variables originales (une variation de log10x pour chaque unité de x est difficile à interpréter).
Bonjour,
oui bien sur c’est possible, le désavantage c’est que parfois la relation n’est pas très simple à expliquer non plus…
Bonjour Claire,
c’est toujours avec un grand plaisir que lis vos articles très pédagogiques et surtout très claires.
Je me posais la question à savoir si les transformation log, sqrt et Box-Cox ne permettaient pas d’avoir un p-value supérieure à 0.05 avec le test de shapiro. Quel(s) modèle(s) devons nous tester.
Merci d’avance de votre réponse
Cordialement
Sidy
Bonjour Sidy,
les transformation log, et Box-Cox, peut être aussi sqrt, permettent d’améliorer la normalité des données, donc d’obtenir, parfois, des pvalues du test de Shapiro-Wilk >0.05.
Bonne continuation
Bonjour Claire,
Merci pour ces précisions
Cordialement
Bonjour,
Article très intéressant et complet, mille merci!
Une petite question: comment choisir entre correlation de Spearman et Kendall?
D’avance, merci pour votre réponse!
Cordialement,
Mathilde
Bonjour Mathilde,
Si mes souvenirs sont bons, je dirais qu’avec Spearman on garde une notion de métrique car c’est basé sur une différence de rangs, alors que Kendall est basé sur les paires concordantes et discordantes.
Les deux méthodes sont applicables. Je pense que Spearman est plus courant que Kendall, mais j’imagine que cela doit aussi dépendre des domaines.
Personnellement, quand la question se pose, je choisis toujours Spearman, mais ce n’est qu’une habitude…
En espérant que cela vous aide.
Bonne continuation.