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

méthode des moindres carrés ordinaires

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") 
Données aberrantes

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") 
Droite de régression en présence de données aberrantes

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") 
Analyse des distance 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") 
Evaluation des outliers par les 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") 
Comparaison des droites de régression avec et sans les points aberrants

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.

régression non paramétrique de theil sen

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.

 

régression quantile

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.

Poursuivez votre lecture

17 réponses

  1. 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

  2. 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

  3. 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 ?

    1. 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 !).

  4. 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

  5. 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à 😉

  6. 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

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée.

Bonjour !

vous venez souvent ?

Identifiez-vous pour avoir accès à toutes les fontionnalités !

Aide mémoire off'R ;)

Enregistrez vous pour recevoir gratuitement mes fiches “aide mémoire” (ou cheat sheets) qui vous permettront de réaliser facilement les principales analyses biostatistiques avec le logiciel R et pour être informés des mises à jour du site.