Régression non paramétrique de Theil Sen
Table des matières
Contexte
La régression non paramétrique de Theil Sen est une alternative à la régression linéaire, qui peut notamment être employée en présence d’outliers.
Pour rappel, la régression linéaire simple est un modèle mathématique, qui met en relation une variable réponse (Y) avec une variable explicative (X), par l’intermédiaire d’une droite.
Cette droite est celle qui passe au plus près de l’ensemble des points. Autrement dit, c’est la droite qui minimise la somme des distances au carré, entre l’observation et son projeté vertical sur la droite. C’est ce qu’on appelle la méthode des moindres carrées (ou MCO pour moindres carrés ordinaires ou encore OLS pour Ordinary Least Squares, en anglais).
Mais parfois les données que l’on étudie comportent quelques observations distantes des autres, qui ne sont pas en accord avec le modèle mathématique sous-jacent. On appelle, parfois ces données, des données aberrantes.
Si ces données s’éloignent suffisamment, alors elles pourront être qualifiées d’outliers, mais ce n’est pas toujours le cas.
# simulation des données
set.seed(1)
n <- 30
x <- round(runif(n=n,min=1, max=10),1)
y <- 2 + 3*x + rnorm(n,0,2.5)
n_outliers <-3
x_out <-c(10,10.5,11)
set.seed(1)
y_out <- 2 + 3*x_out + rnorm(n_outliers,0,2.5) -20
X <- c(x, x_out)
Y <- c(y, y_out)
id <- 1:(n+n_outliers)
type <- c(rep("conforme",n), rep("abberante", n_outliers))
df <- data.frame(id,type, X, Y library(ggplot2)
library(ggrepel)
# visualisation des données
ggplot(df, aes(X, Y, colour=type))+
geom_point(size=3)+
scale_colour_manual(values=c("red", "black"))+
geom_text(aes(label=rownames(df)),size=3,check_overlap = TRUE, vjust=-1) +
theme(legend.position = "bottom") 
Remarque : Les numéros sur le graphique correspondent aux indices (numéros des lignes).
En présence de données aberrantes, la méthode des moindres carrés n’est pas très robuste. Autrement dit, les paramètres de la droite (l’intercept et la pente) sont fortement influencés par ces observations distantes des autres.
Dans cet article, je vous propose d’explorer ce phénomène et quelques approches de modélisation alternatives.
Régression par les moindres carrés
Régression par la méthode des moindre carrés, sur l’ensemble des données
Commençons par ajuster la droite de régression linéaire par la méthode des moindres carrés, sur l’ensemble des observations :
mco <- lm(Y~X, data=df)
coefficients(mco)
## (Intercept) X
## 6.681251 1.936301
confint(mco)
## 2.5 % 97.5 %
## (Intercept) 1.887023 11.475480
## X 1.218787 2.653814 Ici, l’ordonnée à l’origine est égale à 6.68 et la pente à 1.94 avec un intervalle de confiance à 95% compris entre [1.22, 2.65]. Nous pouvons visualiser comment la droite se situe par rapports aux points avec ce code :
ggplot(df, aes(x=X, y=Y))+
geom_point(aes(colour=type))+
scale_colour_manual(values=c("red", "black"))+
geom_text(aes(label=rownames(df), colour=type),size=4,check_overlap = TRUE, vjust=-1)+
geom_smooth(method="lm", se=FALSE, colour="red", size=1)+
annotate("text", x = 3, y = 35, label = "mco : y = 6.7 + 1.9 x ", size=4, colour="red")+
theme(legend.position="bottom") 
Nous voyons que les 3 points aberrants ont attiré la droite de régression vers eux.
Recherche des outliers
On peut alors se demander si ces données aberrantes peuvent être considérées comme des outliers. Si c’est le cas, cela pourrait, éventuellement, dans certaines conditions, constituer un argument pour retirer ces points de l’analyse
Pour cela, nous pouvons employer la fonction check_outliers du package performance qui utilise plusieurs approches de détection des outliers, constitue un score global et fourni une synthèse en sortie.
library(performance)
check_outliers(ols)
## OK: No outliers detected. Ici ces observations ne sont pas considérées comme des outliers. Néanmoins, je regrette un peu que, par défaut, le seuil employé pour les distances de cook soit 1, alors que le seuil généralement employé est 4/n.
Nous pouvons étudier plus précisément les distances de cook et les résidus standardisés de ces 3 points, à l’aide du package broom:
df <- augment(mco)
head(df)
## # A tibble: 6 x 8
## Y X .fitted .resid .hat .sigma .cooksd .std.resid
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 12.1 3.4 13.3 -1.18 0.0560 5.88 0.00130 -0.209
## 2 14.9 4.3 15.0 -0.148 0.0414 5.89 0.0000147 -0.0261
## 3 23.0 6.2 18.7 4.27 0.0304 5.83 0.00881 0.750
## 4 31.7 9.2 24.5 7.16 0.0673 5.73 0.0591 1.28
## 5 11.9 2.8 12.1 -0.218 0.0690 5.89 0.0000565 -0.0390
## 6 31.6 9.1 24.3 7.30 0.0650 5.72 0.0590 1.30
ggplot(df, aes(x=id, y=.cooksd))+
geom_point()+
geom_hline(yintercept=4/nrow(df), linetype="dashed", colour="red")+
ggtitle("Distances de cook") 
Les 3 points ont bien des distances de cook supérieures au seuil. Les distances de cook ne caractérisent pas directement le caractère outlier des données, mais plutôt si celles-ci influencent les résultats de la régression. C’est le cas ici.
Pour le caractère “outliers”, nous pouvons étudier les résidus standardisés. La valeur seuil généralement employée est 3.
ggplot(df, aes(x=id, y=abs(.std.resid)))+
geom_point()+
geom_hline(yintercept=3, linetype="dashed", colour="red")+
ggtitle("Résidus standardisés") 
Nos 3 points sont juste en dessous de ce seuil.
Réalisons, à présent, la régression linéaire par la méthode des moindres carrés, mais en excluant cette fois les 3 points distants.
Remarque : pour obtenir des informations complémentaires sur les outliers, vous pouvez consulter cet article : Comment détecter les outliers avec R .
Régression par les moindres carrés en omettant les 3 points aberrants
Cette analyse va nous permettre de quantifier l’impact de ces 3 points sur les paramètres de la droite de régression :
mco2 <- lm(Y~X, data=df,subset=-c(31,32,33) )
coefficients(mco2)
## (Intercept) X
## 2.279832 2.984457
confint(mco2)
## 2.5 % 97.5 %
## (Intercept) 0.5324377 4.027226
## X 2.7012473 3.267668 L’ordonnée à l’origine passe de 6.68 à 2.28, et la pente de 1.94 à 2.98. Les intervalles de confiance des pentes (avec et sans les 3 points) ne se chevauchent pas, ce qui nous laissent penser que ces estimations sont significativement différentes.
Nous pouvons visualiser ces deux droites, comme ceci :
intercept <- mco2$coefficients[1]
slope <- mco2$coefficients[2]
data_ols2 <- data.frame(X=seq(min(df$X), max(df$X), length.out=100))
data_ols2$Y <- intercept + slope * data_ols2$X
ggplot(df, aes(x=X, y=Y))+
geom_point(aes(colour=type))+
scale_colour_manual(values=c("red", "black"))+
geom_text(aes(label=rownames(df), colour=type),size=4,check_overlap = TRUE, vjust=-1)+
geom_smooth(method="lm", se=FALSE, colour="red", size=1)+
annotate("text", x = 3, y = 35, label = "mco : y = 6.7 + 1.9 x ", size=4, colour="red") +
geom_line(data=data_ols2,aes(x=X, y=Y), colour="blue")+
annotate("text", x = 3, y = 33, label = "mco2 : y = 2.28 + 2.98 x ", size=4, colour="blue")+
theme(legend.position = "bottom") 
Essayons à présent d’autres approches, qui permettent de considérer l’ensemble des observations, mais de rendre les points aberrants moins influants sur la définition de la droite.
Approches alternatives
Parmi ces approches, il ya la régression non paramétrique de Theil Sen, la régression robuste et la régression quantile.
Régression non paramétrique de Theil Sen
La régression non paramétrique de Theil Sen consiste à :
- calculer les pentes des droites entre tous les points deux à deux (il y en a n(n-1)/2, soit ici 528, avec nos 33 points)
- calculer la médiane (m) de ces multiples pentes, il s’agira de l’estimation de la pente de la droite de Theil Sen
- tracer des droites de pente m passant par chacun des points
- calculer la médiane des ordonnées à l’origine de ces multiples droites de pente m ; il s’agira de l’estimation de l’ordonnée à l’origine la pente de la droite de Theil Sen.
La régression non paramétrique de Theil Sen peut être réalisé à l’aide de la fonction mblm(), du package mblm:
library(mblm)
theil_sen<- mblm(Y ~ X,data=df)
coefficients(theil_sen)
## (Intercept) X
## 2.023452 2.962610
confint(theil_sen)
## 0.025 0.975
## (Intercept) 1.959356 4.277086
## X 2.570469 3.005011 Les paramètres de la droite de régression de Theil Sen sont plutôt proches de ceux de la méthode des moindres carrés lorsque les 3 points distants sont retirés.
Plus d’information sur la régression non paramétrique de Theil Sen, ici.
Regresion robuste
Je ne connais pas du tout le principe de cette approche…mais Arthur, en commentaire, nous dit ceci :
“Pour résumer en très grossier (et si j’ai bien compris), on commence par une première recherche des paramètres du modèle en utilisant la méthode des moindres carrées. Une fois le premier modèle fitté, on recommence mais en attribuant un poids à chaque individu (variable X dans ce cas). Le poids associé à un couple (X,Y) est inversement proportionnel à son résidu du modèle précédent. En d’autres termes, si un point est très mal fitté, celui-ci aura moins de poids lors de la prochaine recherche de paramètre par les moindres carrés. Et on recommence le processus x fois jusqu’à temps de trouver une solution qui converge.
C’est évidemment un peu plus compliqué que cela mais l’idée est là”. (Merci Arthur !).
library(MASS)
#fit robust regression model
robust <- rlm(Y~X, data=df)
summary(robust)
##
## Call: rlm(formula = Y ~ X, data = df)
## Residuals:
## Min 1Q Median 3Q Max
## -20.69260 -1.54148 -0.02664 1.69417 3.46758
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 3.2547 1.0335 3.1492
## X 2.7590 0.1547 17.8375
##
## Residual standard error: 2.407 on 31 degrees of freedom
coefficients(robust)
## (Intercept) X
## 3.254700 2.758985 Là encore, la droite obtenue est relativement proche de la droite des moindres carrés sans les 3 points distants et de la droite de Theil Sen.
Quantile regression
La méthode des moindres carrés fournit une estimation de la moyenne de la variable réponse (Y) en fonction d’une valeur de la variable explicative (X) prédictives.
La régression quantile, quant à elle, donne une estimation d’un quantile donné de la réponse (la médiane si le quantile =0.5) en fonction d’une valeur de la variable explicative (X).
La pente estimée représente l’évolution du quantile concerné de la variable réponse lorsque la variable explicative augmente d’une unité.
Nous allons essayer en choisissant le quantile 0.5, ce qui devrait nous permettre d’obtenir une droite assez proche de celle de Theil Sen
library(quantreg)
quantreg <- rq(Y ~ X, data = df, tau=.5)
summary(quantreg)
##
## Call: rq(formula = Y ~ X, tau = 0.5, data = df)
##
## tau: [1] 0.5
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 1.92259 1.24060 4.40125
## X 2.98973 2.07993 3.18376
coefficients(quantreg)
## (Intercept) X
## 1.922588 2.989729 C’est effectivement le cas.
Synthèse
Les 3 méthodes alternatives de la régression linéaire conduisent à des droites de régression plutôt similaires et très proche de la droite obtenue par la méthode des moindres carrés en omettant les points aberrants.
Parmi toutes ces méthodes alternatives, la régression de Theil Sen me parait la plus intéressante, car elle est facile à comprendre, et assez intuitive.
Et vous, qu’en pensez-vous ? Dites-le moi en commentaires
Si cet article vous a plu, ou vous a été utile, vous pouvez soutenir le blog, en réalisant un don libre sur sa page Tipeee.
18 Responses
Bon travail Claire ! Toujours super documenté et illustré.
A noter les packages robust et robustbase que j’apprécie.
Et plus encore dans https://cran.r-project.org/web/views/Robust.html
Bonjour Samuel,
Merci également pour les références de ces packages !
Bonjour Claire,
Merci beaucoup pour cet article, bien expliqué. La régression de Theil-Sen est en effet un bon moyen de tenir compte des outliers.
À ce sujet, il y a un package appelé Openair (développé par David Carslaw) basé sur la méthode Theil-Sen, qui tient également compte de l’autocorrélation du signal dans le calcul des intervalles de confiance. Est-ce qu’on pourrait avoir un petit article explicatif là-dessus de votre part, si vous avez un peu de temps à y consacrer ? =)
Yann
Bonjour Yann,
merci pour la référence. Je note la suggestion d’article.
Bien à vous.
cet article tombe pile au moment où je suis en discussion avec un de mes encadrants sur ce que je dois faire des outliers présents dans mon jeu de données.
merci Claire
Bonjour CLAIRE, Très ravis de cette publication. Est-il possible de faire la régression linéaire multiple dans de cas non paramétrique. Càd corréler une variable quantitative avec autant des variables explicatives.
Merci
Tres bonne approche.
Vous etes a remercier infiniment avec toutes ces fouilles.
Bonsoir Claire, merci pour cette découverte, je ne connaissais cette régression non paramétrique. Question : ls moindres carrés ne valent que si on fait l’hypothèse d’un modèle d’erreur Gaussien sur Y. Une régression non paramétrique permet a priori de s’affranchir de cette hypothèse. Du coup, peut-on l’utiliser quelle que soit la nature des données Y ? Binaires, comptage, quantitative continues ?
Bonjour Sandrine,
D’après ma compréhension, pour estimer les paramètres de la droite (pente et intercept), la méthode des moindres carrés ne fait aucune hypothèse sur la distribution des erreurs. Par contre, c’est l’estimation de l’erreur standard des paramètres, qui découle de la méthode des moindres carrés qui repose sur une distribution Gaussienne des erreurs. Dans les sorties de l’utilisation de la méthode de Theil Sen, on obtient un intervalle de confiance. Je ne sais pas comment il est calculé, mais sans doute par bootstrap. Du coup, je pense qu’on peut utiliser cette approche quelle que soit la nature des données (je ferais quand même quelques petits tests pour vérifier !).
Super travail et article pédagogique (comme souvent…. NON : comme toujours ! 🙂 ). A bientôt.
Bonjour ,je me jouis de cet immense travail.toutes mes félicitations à vous.
Quel logiciel utilisez-vous pour ces calculs ?
Bonjour,
vraiment, vous n’avez pas une petite idée du logiciel ??
C’est vraiment un article intéressant et innovant pour nous les débutants en R. C’était toujours confronté à des tels cas mais je ne savais m’y prendre
Merci
Cordialement
Bonjour Claire,
Merci beaucoup pour cet article très intéressant et qui tombe à pic !
Concernant la régression robuste, la fonction rlm utilise la méthode “iterated re-weighted least squares” pour fitter les coefficients du modèle (https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares).
Pour résumer en très grossier (et si j’ai bien compris), on commence par une première recherche des paramètres du modèle en utilisant la méthode des moindres carrées. Une fois le premier modèle fitté, on recommence mais en attribuant un poids à chaque individu (variable X dans ce cas). Le poids associé à un couple (X,Y) est inversement proportionnel à son résidu du modèle précédent. En d’autres termes, si un point est très mal fitté, celui-ci aura moins de poids lors de la prochaine recherche de paramètre par les moindres carrés. Et on recommence le processus x fois jusqu’à temps de trouver une solution qui converge.
C’est évidemment un peu plus compliqué que cela mais l’idée est là 😉
Bonjour Arthur,
Merci pour l’explication. Je vais l’ajouter à l’article.
merci bcp Claire, une très bonne alternative pour les cas difficiles. avant cet article, je savait juste qu’elle existe, mais avec votre article , je vois plus clair et avec l’outil R c’est juste génial!
merci et bon courage
Article très intéressant et utile ! J’adore apprendre de nouvelles méthodes en stat…
Merci beaucoup pour tout ce que vous faites pour rendre ces approches bien claires pour moi. Les données quantitatives contiennent souvent des valeurs abberantes qu’il faut savoir les traiter convenablement pour avoir des estimations moins biaisées. J’en tiendrai compte pour mes travaux de modélisation.