Comment détecter les outliers avec R

comment lire un boxplot

Table des matières

C'est quoi un outlier ?

Un outlier, ou donnée aberrante est « une valeur ou une observation qui est « distante » des autres observations effectuées sur le même phénomène, c’est-à-dire qu’elle contraste grandement avec les valeurs « normalement » mesurées. Une donnée aberrante peut être due à la variabilité inhérente au phénomène observé ou bien elle peut aussi indiquer une erreur expérimentale. Les dernières sont parfois exclues de la série de données ».(wikipedia)

Détection des outliers par boxplots ?

Une méthode classiquement employée pour détecter les outliers, consiste à réaliser un boxplot. On parle alors de méthode de détection univariée car elle ne concerne qu’une seule dimension, ou variable.

Réalisation des boxplots

Voici comment réaliser des boxplots avec le package ggplot2. Je prends ici comme exemple les données « mpg » contenues dans le package « ggplot2 », et je m’intéresse à la variable « hwy » (Highway miles per gallon) en fonction du type de véhicule (variable class)

library(ggplot2)
ggplot(mpg, aes( x=class,y=hwy, fill=class)) +
    geom_boxplot()+ 
    xlab(label = "Type of car") +
    ylab(label = "Highway miles per gallon") +
    theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1))+
    theme(legend.position="none")+
    ggtitle("Exemple de boxplots sur les données mpg") 
outlier avec R

Sur cette visualisation des données, les outliers sont représentés sous forme de points. Ils correspondent à des observations dont les valeurs sont :

  • supérieures à la valeur du 3ème quartile plus 1.5 fois l’intervalle inter-quartile,
  • ou inférieures à la valeur du 1er quartile moins 1.5 fois l’intervalle inter-quartile.

Règle d'interprétation des boxplots

Pour celles et ceux qui ne sont pas très familiers avec les boxplots, voici un rappel des règles d’interprétation (d’après Laura DeCicco).

interprétation boxplot

Récupération des valeurs outliers

Imaginons, par exemple, que je souhaite obtenir les valeurs des observations outliers (concernant la variable hwy), correspondant aux voitures de type SUV (le boxplot rose sur grap du paragraphe 2.1). Pour cela, j’utilise la fonction boxplot.stats() :

# sélection des données suv
suv <- mpg %>%
    filter(class=="suv")
boxplot.stats(suv$hwy)$out

## [1] 12 12 25 24 27 25 26 23 

Récupérer l'indice des outliers

Il peut être intéressant de connaître sur quelles lignes du jeu de données se trouvent les outliers. Pour cela, on peut garder les valeurs outliers dans un objet, puis l’utiliser en argument de la fonction which(), comme ceci :

outlier_val <- boxplot.stats(suv$hwy)$out
outlier_val

## [1] 12 12 25 24 27 25 26 23
outlier_idx <-which(suv$hwy %in% c(outlier_val))
outlier_idx

## [1] 12 30 49 50 51 52 53 54 

Pour vérifier, on peut afficher les lignes correspondant à outlier_idx:


suv[outlier_idx,]

    ## # A tibble: 8 x 11
    ##   manufacturer model displ  year   cyl trans drv     cty   hwy fl    class
    ##                    
    ## 1 dodge        dura~   4.7  2008     8 auto~ 4         9    12 e     suv  
    ## 2 jeep         gran~   4.7  2008     8 auto~ 4         9    12 e     suv  
    ## 3 subaru       fore~   2.5  1999     4 manu~ 4        18    25 r     suv  
    ## 4 subaru       fore~   2.5  1999     4 auto~ 4        18    24 r     suv  
    ## 5 subaru       fore~   2.5  2008     4 manu~ 4        20    27 r     suv  
    ## 6 subaru       fore~   2.5  2008     4 manu~ 4        19    25 p     suv  
    ## 7 subaru       fore~   2.5  2008     4 auto~ 4        20    26 r     suv  
    ## 8 subaru       fore~   2.5  2008     4 auto~ 4        18    23 p     suv 

Juxtaposer les observations sur les boxplots

Lorsque je réalise un boxplot, j’ai pour habitude de faire également figurer les observations, parce que ça permet de se faire une meilleure idée du nombre de données et de leur répartition. Pour cela, j’utilise la couche geom_jitter().

Remarque : L’argument width de la fonction geom_jitter() permet de gérer la largeur de répartition des points à l’intérieur du boxplot.

ggplot(mpg, aes( x=class,y=hwy, fill=class, colour=class)) +
    geom_jitter(width=0.25) +
    geom_boxplot(alpha=0.5)+ 
    xlab(label = "Type of car") + 
    ylab(label = "Highway miles per gallon") + 
    theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1))+
    theme(legend.position="none")+
    theme_classic()+ 
    ggtitle("Boxplot avec les observations") 
boxplot outlier avec R

Le problème de ce type de boxplot, c’est que les outliers sont représentés deux fois: une fois en tant qu’observation et une fois en tant qu’outlier (les points sont plus clairs, ils sont particulièrement visibles dans la partie supérieure du boxplot rose).

Pour éviter cette double représentation, j’utilise l’argument outlier.shape=NA, ce qui permet de ne pas afficher la version outlier :


 ggplot(mpg, aes( x=class,y=hwy, fill=class, colour=class)) +
        geom_jitter(width=0.25)+
        geom_boxplot(alpha=0.5, outlier.shape=NA)+  
        xlab(label = "Type of car") +
        ylab(label = "Highway miles per gallon") +
        theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1))+
        theme(legend.position="none")+
        theme_classic()+
        ggtitle("Boxplot avec les observations") 
boxplot outlier avec ggplot2

A présent, les outliers sont simplement tous les points au delà des traits verticaux de part et d’autre de chaque boite.

Les autres méthodes de détection univariées

Méthode basée sur les percentiles

Il existe d’autres méthodes de détection des outliers. L’une d’elles consiste à considérer comme outlier les données se situant à l’extérieur de l’intervalle constituté par les percentiles 2.5 et 97.5 par exemple, ou encore par les percentiles 1 et 99.

pinf = 0.025 # fixe le percentile de limite inférieure 
psup = 0.97 5# fixe le percentile de limite inférieure 
binf <- quantile(suv$hwy,pinf) # calcule la borne inf de l'intervalle 
binf 
## 2.5% 
## 13.05 
bsup <- quantile(suv$hwy,psup) # calcule la borne sup de l'intervalle bsup 
## 97.5% 
## 25.475 

On peut obtenir les indices des lignes des outliers comme ceci:

outlier_idx <- which(suv$hwy < binf | suv$hwy > bsup) 
outlier_idx 
## [1] 12 30 51 53 

Puis leur valeur :

outlier_val <- suv[outlier_idx,"hwy"] 
outlier_val 
## # A tibble: 4 x 1 
## hwy 
## 
## 1 12 
## 2 12 
## 3 27 
## 4 26 

La méthode de Hampel

Une autre méthode, dite de Hampel, consiste à considérer comme outliers les valeurs en dehors de l’intervalle constitué par la médiane, plus ou moins 3 déviation absolue de médiane :

\[ I=[median – 3\;\ast mad \;;\;median + 3\;\ast mad ]  \]

Avec : 

\[ mad = \text{Median Absolute Deviation} \]

Et:

\[  mad=median\;(|y_{i} – \tilde{y} |) \]

et

\[ \tilde{y} = median(y) \]

Parfois, une autre valeur que 3 est employée pour calculer l’intervalle.

Voici comment constituer l’intervalle selon la méthode de Hampel.

k= 3 
# calcule la borne inf de l'intervalle binf
binf <- median(suv$hwy) - k * mad (suv$hwy) 
binf 
## [1] 10.8283 

# calcule la borne sup de l'intervalle bsup 
bsup <- median(suv$hwy) + k * mad (suv$hwy) 
bsup
## [1] 24.1717 

Il est ensuite facile d’obtenir l’indice des lignes contenant les outliers :.

outlier_idx <- which(suv$hwy < binf | suv$hwy > bsup)
outlier_idx

## [1] 49 51 52 53

Puis leur valeur :

outlier_val <- suv[outlier_idx,"hwy"]
outlier_val

## # A tibble: 4 x 1
## hwy
## 
## 1 25
## 2 27
## 3 25
## 4 26 

Le test statistique de Grubbs

Le test statistique de Grubbs permet de tester si la valeur la plus élevée est un outlier. Son hypothèse nulle spécifie que la valeur la plus élevée n’est pas un outlier, alors que son hypothèse alternative spécifie que la valeur la plus élevée est un outlier. Si la p-value du test est inférieure au seuil de significativité choisi (en général 0.05) alors on concluera que la valeur la plus élevée est outlier. Pour réaliser ce test avec R, on utilise la fonction grubbs.test() du package « outliers » :

library(outliers) 
grubbs.test(suv$hwy)
## 
## Grubbs test for one outlier 
## 
## data: suv$hwy 
## G = 2.97890, U = 0.85215, p-value = 0.06296 
## alternative hypothesis: highest value 27 is an outlier 

Ici la p-value est proche du seuil de significativité, j’aurais tendance à considérer la valeur la plus haute (27) comme étant un outlier.

Il est possible de tester le caractère outlier de la plus faible valeur en utilisant l’argument opposite=TRUE :

grubbs.test(suv$hwy, opposite = TRUE) 
## 
## Grubbs test for one outlier 
## 
## data: suv$hwy 
## G = 2.05810, U = 0.92942, p-value = 1 
## alternative hypothesis: lowest value 12 is an outlie 

Ici, la plus faible valeur (12) n’est pas considérée comme un outlier.

Détection des outliers sur les résidus

Dans le cadre d’une régression, il est possible de mettre en évidence des outliers sur les résidus, en employant la fonction outlierTest() du package « car ». Je réalise ici une régression linéaire simple de la variable « hwy » (highway miles per gallon) en fonction de la variable « displ » (engine displacement, in litres) :

ggplot(suv, aes(x=displ, y=hwy))+
   geom_point()+
   geom_smooth(method="lm") 
outliers sur résidus
fit <- lm(hwy~displ, data=suv) 

Les plots des distances de cook  permettent déjà d’identifier certains points outliers, ici les données 51, 53, et 49.

plot(fit,4) 
distance de cook

library(car)
outlierTest(fit)

## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferonni p
## 12 -2.644498 0.010466 0.64891 

Ici, la p-value ajustée par la méthode de Bonferonni est très éloignée du seuil de 0.05. Le résidus le plus élevé, ne peut donc pas être considéré comme outlier. Pour vérifier quel est l’indice du résidu le plus élevé, on peut utiliser la commande suivante :

which.max(rstudent(fit))

## 51
## 51 

Que faire des outliers ?

En premier lieu, les données détectées comme outliers doivent être vérifiées.

Pour la suite, il n’y a pas de réponse unique à cette question !

Dans certains domaines, comme en génomique ou protéomique, il n’est pas rare de supprimer les outliers du jeu de données, car ces données extrêmes sont considérées comme issues d’un processus qui a dysfonctionné.

Dans d’autres domaines, au contraire, ces données extrêmes sont conservées, mais les analyses statistiques sont conduites avec et sans, afin d’évaluer leur impact sur les résultats, voir un exemple dans le paragraphe 7 de cet article .

Pour aller plus loin

Vous trouverez d’autres méthodes d’identification des outliers dans le chapitre 7 du livre « R and Data Mining: Examples and Case Studies« .

Le package « Outliers » contient également plusieurs fonctions dédiées à la détection et à la gestion des outliers.

ls("package:outliers") 
## [1] "chisq.out.test" "cochran.test" "dixon.test" "grubbs.test" 
## [5] "outlier" "pcochran" "pdixon" "pgrubbs" 
## [9] "qcochran" "qdixon" "qgrubbs" "qtable" 
## [13] "rm.outlier" "scores" 

Le code des règles de décision du boxplot

J’ai adapté le code original développé par Laura DeCicco et disponible là.

library(tidyverse) 
library(cowplot) 
ggplot_box_legend <- function(family = "serif"){ 
# Create data to use in the boxplot legend: set.seed(100)
   sample_df <- data.frame(parameter = "test", values = sample(500)) 
# Extend the top whisker a bit:
sample_df$values[1:100] <- 701:800 
# Make sure there's only 1 lower outlier: 
 sample_df$values[1] <- -350 
# Function to calculate important values: 
  ggplot2_boxplot <- function(x){ quartiles <- as.numeric(quantile(x, probs = c(0.25, 0.5, 0.75))) 
  names(quartiles) <- c("25ième percentile", "50ième percentilen(médiane)", "75ième percentile") 
  IQR <- diff(quartiles[c(1,3)]) 
  upper_whisker <- max(x[x < (quartiles[3] + 1.5 * IQR)]) 
  lower_whisker <- min(x[x > (quartiles[1] - 1.5 * IQR)]) 
  upper_dots <- x[x > (quartiles[3] + 1.5*IQR)] 
  lower_dots <- x[x < (quartiles[1] - 1.5*IQR)] 
  return(list("quartiles" = quartiles, "25ième percentile" = as.numeric(quartiles[1]), "50ième percentilen(médiane)" = as.numeric(quartiles[2]), "75ième percentile" = as.numeric(quartiles[3]), "IQR" = IQR, "upper_whisker" = upper_whisker, "lower_whisker" = lower_whisker, "upper_dots" = upper_dots, "lower_dots" = lower_dots)) } 
  # Get those values:
    ggplot_output <- ggplot2_boxplot(sample_df$values) 
    # Lots of text in the legend, make it smaller and consistent font: 
    update_geom_defaults("text", list(size = 3, hjust = 0, family = family)) 
    # Labels don't inherit text:
    update_geom_defaults("label", list(size = 3, hjust = 0, family = family)) 
    # Create the legend: 
    # The main elements of the plot (the boxplot, error bars, and count) are the easy part. 
    # The text describing each of those takes a lot of fiddling to 
    # get the location and style just right: 
    explain_plot <- ggplot() + stat_boxplot(data = sample_df,aes(x = parameter, y=values),geom ='errorbar', width = 0.3) + 
    geom_boxplot(data = sample_df,aes(x = parameter, y=values), width = 0.3, fill = "orange") + 
    # geom_text(aes(x = 1, y = 950, label = "500"), hjust = 0.5) + 
    # geom_text(aes(x = 1.17, y = 950, label = "Nombre d'observations"),fontface = "bold", vjust = 0.4) +
    theme_minimal(base_size = 5, base_family = family) +
    geom_segment(aes(x = 2.3, xend = 2.3, y = ggplot_output[["25ième percentile"]], yend = ggplot_output[["75ième percentile"]])) +
    geom_segment(aes(x = 1.2, xend = 2.3, y = ggplot_output[["25ième percentile"]], yend = ggplot_output[["25ième percentile"]])) +
    geom_segment(aes(x = 1.2, xend = 2.3, y = ggplot_output[["75ième percentile"]], yend = ggplot_output[["75ième percentile"]])) +
    geom_text(aes(x = 2.4, y = ggplot_output[["50ième percentilen(médiane)"]]), label = "Intervallen interquartile", fontface = "bold",vjust = 0.4) +
    geom_text(aes(x = c(1.17,1.17), y = c(ggplot_output[["upper_whisker"]], ggplot_output[["lower_whisker"]]), label = c("Plus grande valeur dans l'étendue 1.5 foisn l'intervalle interquartile au dessusn du 75ième percentile", "Plus petite valeur dans l'étendue 1.5 foisn l'intervalle interquartile au dessousn du 25ième percentile")), fontface = "bold", vjust = 0.9) + 
    geom_text(aes(x = c(1.17), y = ggplot_output[["lower_dots"]], label = "Valeur Outlier "), vjust = 0.5, fontface = "bold") + geom_text(aes(x = c(1.9), y = ggplot_output[["lower_dots"]],label = " La valeur est >1.5 fois et"), vjust = 0.5) + geom_text(aes(x = 1.17, y = ggplot_output[["lower_dots"]], label = "<3 fois l'intervalle interquartilen au delà de chaque coté de la boite"), vjust = 1.5) +
    geom_label(aes(x = 1.17, y = ggplot_output[["quartiles"]], label = names(ggplot_output[["quartiles"]])), vjust = c(0.4,0.85,0.4), fill = "white", label.size = 0) + 
    ylab("") +
    xlab("") + 
    theme(axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), aspect.ratio = 4/3, plot.title = element_text(hjust = 0.5, size = 10)) +
    coord_cartesian(xlim = c(1.4,3.1), ylim = c(-600, 900)) + 
    labs(title = "INTERPRETATION") return(explain_plot) } 

Conclusion

Et vous, quelle méthode utilisez vous pour détecter les outliers ? Et comment les gérer vous dans la suite de vos analyses statistiques ?

Et comme d’habitude, si cet article vous a plus, ou vous a été utile, et si vous le souhaitez, vous pouvez soutenir ce blog en faisant un don sur sa page Tipeee 

Credits photos : geralt

10 Responses

  1. Bonjour Claire, très intéressant l’article, il va nous être très utile car c’est notre prochain objectif : valider les données !
    Merci et à bientôt
    Miriam BASSO

  2. Bonjour,
    je voudrais vous souhaiter une bonne année.
    je suis entrain de préparer une publication sur l’écologie des dinoflagellés épiphytes. mes résultats comportent des outliers. votre article va m’aider beaucoup dans l’interprétation de mes graphes.

    merci

  3. Bonjour ,
    j’ai aimé votre publication. je voudrais juste apporter juste une petite extension du test de grubbs. En fait ce test est itératif. par exemple si la plus grande valeur est aberrante, il faut continuer à test si la valeur qui vient avant cette valeur est aussi aberrante. vous pouvez revoir ce code.

    compact% filter(class== »compact »)
    nrow(compact)
    grubbs.test(compact$hwy) #la born sup est abberant
    ##test de grubbs itératif
    #names(grubbs.test(compact$hwy))
    compact_bis<-compact
    while(grubbs.test(compact$hwy)$p.value <0.05){
    compact% filter(hwy<max(hwy,na.rm = T))
    }
    nrow(compact)
    nrow(compact_bis)
    outlier2 <-setdiff(compact_bis$hwy,compact$hwy) #"" compar 2 vecteu
    compact_bis% mutate(out=hwy %in% outlier2)
    compact_bis

  4. Bonjour,

    Merci beaucoup pour votre article !!

    Petite question complémentaire, vous avez choisi de connaître les outliers uniquement de la catégorie SUV de votre variable class.
    Si on souhaite obtenir les outliers de chaque catégorie de cette variable, quelle formule peut-on utiliser ? (A part filtrer une à une les catégories…)

    En vous remerciant par avance,

  5. Bonjour Claire,
    Comme toujours, vos articles sont très qualitatifs et apportent beaucoup d’informations pertinentes. Merci à vous pour ces publications.
    En ce qui concerne la détection de valeurs aberrantes, je suis amené à les recherches dans le cadre des validations des méthodes d’analyse. Dans ce contexte, il semblerait que les tests privilégiés soient principalement les tests de Cochran et de Grubbs.
    En ce qui concerne le test de Grubbs, on peut en mettre en oeuvre deux types : un pour rechercher une valeur extrême seulement (test de Grubbs simple) ou une paire de valeurs extrêmes (test de Grubbs double au cas où la valeur extrême serait « camouflée » par une valeur trop proche d’elle). Il y a aussi un test de Grubbs double qui teste les 2 valeurs extrêmes (la plus faible et la plus grande) à la fois (mais j’avoue ne pas bien comprendre l’intérêt de ce dernier).

    De plus, comme les valeurs extrêmes sont toujours « délicates » à traiter (une donnée obtenue devant essentiellement être considérée comme existante et comprendre donc pourquoi elle est apparue…), il convient d’être prudent avant de la considérer comme aberrante. Ainsi, dans le protocole proposé pour les validations des methodes d’analyse, on considéra qu’une valeur est aberrante à partir du moment où la p-value du test renvoie une valeur inférieure à 1%. Si cette p-value est seulement inférieure à 5% (et supérieure à 1%), on considère alors que la valeur testée est suspecte.
    La décision finale quant au devenir de cette valeur extrême reste néanmoins floue ou du moins à la discrétion de l’expérimentateur (on la supprime en général lorsqu’on la détectée comme aberrante mais la décision est laissée au choix lorsqu’elle n’est que suspecte)…

    Bien cordialement

    1. Bonjour Philippe,

      Merci pour votre commentaire. Depuis la rédaction de cet article, j’ai découvert le package performance qui propose une fonction check_outliers(), qui permet d’utiliser plusieurs méthodes (celle de Tukey avec les IQR, celle bassée sur les percentiles etc..), il y a même des méthodes mutivariées.

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée. Les champs obligatoires sont indiqués avec *

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.