Traitement des données d’enquêt ENSAE

Chapitre 1: Traitement des données manquantes

Introduction

Il arrive assez fréquemment que à l’issu d’une enquête des observations soient incomplètes, c’est à dire les valeurs d’une ou plusieurs variables manquent. On parle d’observations à données manquantes. Par exemple, dans le cas d’un sondage, certaines personnes ne répondent pas à certaines questions ou alors certaines réponses sont oubliées lors de la saisie informatique. Quel est l’impact de l’absence de certaines données sur la construction de modèles à partir des données et sur la prise de décision avec ces modèles ?

Cause des données manquantes

  • un refus total ou partiel (surtout pour les questions jugées indiscrètes)

  • une omission ou une erreur de saisie;

  • Impossibilité de contacter l’unité ou une perte de vue.

  • Perte de questionnaire ou questionnaire inexploitable.

  • Abandon de l’enquêté (pour les questionnaires longs ou les enquêtes WEB) Données aberrantes;

Conséquences des données manquante

  • Introduction d’un biais de non réponse (Ex: sous-estimation: revenu, viols etc.)
  • Perte de précisions;
  • Diminution de la taille de l’échantillon.

Caractéristique des données manquantes

Afin de mieux comprendre l’impact de différentes approches proposées pour faire face au problème des données manquantes, il est nécessaire d’examiner d’abord la nature de ce manque de données.

  • Données manquant de façon complètement aléatoire (missing completely at random, MCAR)
  1. La probabilité qu’une variable 𝑋_𝑖 soit manquant est indépendante des valeurs prises les autres variables.
  2. On ne peut pas établir dans ce cas le profil des individus ayant une valeur manquante
  • Données manquant de façon aléatoire (missing at random, MAR)
  1. La probabilité de la DM peut dépendre des valeurs des autres variables.
  • Données manquant de façon non aléatoire (missing not at random, MNAR)
  1. donnée manquante est du aux valeurs non observées des autres variables.

Correction des données manquantes

Méthodes d’imputations

  • Compléter la non réponse par une valeur plausible.
  • L’imputation des données manquantes : pas très simple
  • Existence de plusieurs méthodes;
  • Traiter cas par cas (Pas méthodes susceptibles de détecter toutes les incohérences à priori, Ou de prévoir à l’avance les méthodes d’imputations les plus pertinentes).

Methodes

On note différentes méthodes d’imputation:

  • méthode déductive : on déduit la réponse à partir des valeurs observées pour les autres variables.

  • Cold-deck : on remplace la DM par une information provenant d’une autre source de données.

  • Imputation par la moyenne ou la moyenne par classe : On remplace les valeurs manquantes par la moyenne observée chez les répondants ou un groupe de répondant. [on peut prendre aussi la médiane à la place]

  • Imputation par régression: Modèle de régression dépendra des cas (logit, GLM, AFD, Tobit, etc.)

  • Imputation par tendance unitaire: Pour les enquêtes répétées dans le temps (présence de série chronologiques). On multiplie la valeur précédente par un facteur d’actualisation qui peut être calculé sur les répondants pour les répondants de la même classe.

  • Hot-deck simple ou hot-deck par classe : On remplace la valeur manquante par une valeur choisi aléatoirement dans l’échantillon ou dans un groupe de l’échantillon;

  • Hot-deck séquentielle ou par classe: On remplace la valeur manquante par une valeur de l’individu précédent direct ou ayant des caractéristiques similaires ;

  • Hot-deck métrique: On remplace par la valeur de l’individu le plus proche à l’aide d’une mesure de distance (plus proche voisin) ;

  • Imputation par les méthodes d’analyses factorielles (ACP, AFC et ACM)

Cas pratique

Traitement de données manquantes partielles

library(haven)
library(dplyr)
## 
## Attachement du package : 'dplyr'
## Les objets suivants sont masqués depuis 'package:stats':
## 
##     filter, lag
## Les objets suivants sont masqués depuis 'package:base':
## 
##     intersect, setdiff, setequal, union
library(questionr)
## Warning: le package 'questionr' a été compilé avec la version R 4.2.3
library(readxl)
producteurs <- read_excel("D:/OneDrive - CETUD/Bureau/cours OD/ENSAE/Donnees enquete marches betail.xlsx")
glimpse(producteurs)
## Rows: 332
## Columns: 9
## $ Espèces               <chr> "BOVINS", "BOVINS", "BOVINS", "BOVINS", "BOVINS"…
## $ Catégories            <chr> "Taureau", "Bœuf", "Taurillon", "Vache", "Géniss…
## $ `Effectifs présentés` <dbl> 78, 88, 428, 503, 209, 98, 88, 458, 541, 208, 52…
## $ `Effectifs vendus`    <dbl> 59, 73, 345, 279, 186, 45, 37, 302, 322, 106, 40…
## $ `Effectifs exportés`  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 101, 95,…
## $ `Prix moyens`         <dbl> 458333, 400050, 128333, 106667, 188333, 60000, 5…
## $ Marché                <chr> "Altou Pass", "Altou Pass", "Altou Pass", "Altou…
## $ Région                <chr> "Tambacounda", "Tambacounda", "Tambacounda", "Ta…
## $ Année                 <dbl> 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, …
describe(producteurs)
## [332 obs. x 9 variables] tbl_df tbl data.frame
## 
## $Espèces: 
## character: "BOVINS" "BOVINS" "BOVINS" "BOVINS" "BOVINS" "BOVINS" "BOVINS" "OVINS" "OVINS" "OVINS" ...
## NAs: 0 (0%) - 7 unique values
## 
## $Catégories: 
## character: "Taureau" "Bœuf" "Taurillon" "Vache" "Génisse" "Veau" "Velle" "Mâle" "Femelle" "Mâle castré" ...
## NAs: 0 (0%) - 23 unique values
## 
## $Effectifs présentés: 
## numeric: 78 88 428 503 209 98 88 458 541 208 ...
## min: 1 - max: 4150 - NAs: 7 (2.1%) - 156 unique values
## 
## $Effectifs vendus: 
## numeric: 59 73 345 279 186 45 37 302 322 106 ...
## min: 0 - max: 2100 - NAs: 16 (4.8%) - 134 unique values
## 
## $Effectifs exportés: 
## numeric: NA NA NA NA NA NA NA NA NA NA ...
## min: 1 - max: 300 - NAs: 257 (77.4%) - 43 unique values
## 
## $Prix moyens: 
## numeric: 458333 400050 128333 106667 188333 60000 56667 118333 46667 54167 ...
## min: 650 - max: 922500 - NAs: 0 (0%) - 265 unique values
## 
## $Marché: 
## character: "Altou Pass" "Altou Pass" "Altou Pass" "Altou Pass" "Altou Pass" "Altou Pass" "Altou Pass" "Altou Pass" "Altou Pass" "Altou Pass" ...
## NAs: 0 (0%) - 23 unique values
## 
## $Région: 
## character: "Tambacounda" "Tambacounda" "Tambacounda" "Tambacounda" "Tambacounda" "Tambacounda" "Tambacounda" "Tambacounda" "Tambacounda" "Tambacounda" ...
## NAs: 0 (0%) - 11 unique values
## 
## $Année: 
## numeric: 2019 2019 2019 2019 2019 2019 2019 2019 2019 2019 ...
## min: 2019 - max: 2019 - NAs: 0 (0%) - 1 unique values
library(VIM)
## Warning: le package 'VIM' a été compilé avec la version R 4.2.3
## Le chargement a nécessité le package : colorspace
## Le chargement a nécessité le package : grid
## VIM is ready to use.
## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues
## 
## Attachement du package : 'VIM'
## L'objet suivant est masqué depuis 'package:datasets':
## 
##     sleep
library(mice)
## Warning: le package 'mice' a été compilé avec la version R 4.2.3
## 
## Attachement du package : 'mice'
## L'objet suivant est masqué depuis 'package:stats':
## 
##     filter
## Les objets suivants sont masqués depuis 'package:base':
## 
##     cbind, rbind

imputation par moyenne simple

Ici on cherche à remplacer les valeurs manquantes sur la variable Effectifs exportés par la moyenne des effectifs exportés

producteurs$eff_expor<-producteurs$'Effectifs exportés'

 #On s'intéresse à la variable effectif exportés 
producteurs<-producteurs %>% mutate(eff_expor_moy=if_else(is.na(eff_expor), floor(mean(eff_expor, na.rm=T)), floor(eff_expor)))
View(producteurs)

Imputation par la moyenne par région

Ici on également à imputer la valeur manquante par la moyennes des effectifs exportés par région

producteurs$region<-producteurs$Région
producteurs<-producteurs %>% group_by(region) %>% mutate(eff_expor_moy_reg=if_else(is.na(eff_expor), mean(eff_expor, na.rm=T),eff_expor))

Hotdeck simple

producteurs$eff_expor_hd <- producteurs$eff_expor
producteurs <- hotdeck(producteurs, variable = 'eff_expor_hd')

Hotdeck par domaine

producteurs$eff_expor_hdd <- producteurs$eff_expor

producteurs <- hotdeck(producteurs, variable = 'eff_expor_hdd', domain_var ="region")

Hotdeck séquentiel

## producteurs <- hotdeck(producteurs, variable = 'prod_hdd', domain_var =c("regions", "ja13"), ord_var = c("superficie", "rendement"))

KNN

##producteurs$prod_kNN<- producteurs$production_kg

##producteurs <- kNN(producteurs, variable = 'prod_kNN', dist_var = c("superficie","ja16", "ja18"), k=5 )

Chapitre 2: Traitement des données abérantes

Qu’est-ce qu’une valeur aberrante ou atypique?

  • Une observation ou ensemble d’observations qui diffère de façon significative de la tendance globale des autres observations.

  • Une observation incohérente par rapport aux autres données de l’ensemble ou qui ne semble pas respecter une norme ou une relation bien définie.

Elle peut être unidimensionnelle ou multidimensionnelle

Conséquence

  • Augmentation de la variance des estimateurs;

  • Estimateurs biaisés (moyenne, étendu, etc.)

  • Déformation des nuages de variables et individus en analyse bi-variée ou factorielle

Cause

  • Caractère aléatoire des sondages: une variabilité inhérente dans la population;

  • Erreur de mesure: erreur de collecte ou de saisie

Détection :

  • Contrôle uni-varié : se fixer une borne maximale ou minimale.

  • Détection graphique : box plot, nuage de points, histogramme etc.

  • Contrôle logique: nuage de point ou tableau croisé avec une variable corrélée.

  • Détermination de seuil : On cherche dans les valeurs en dehors de l’intervalle \[[\bar{X}−1,5(𝑄_3−𝑄_1 );\bar{X}+1,5(𝑄_3−𝑄_1)]\]

  • Analyse multivariée

Détection- Test de grubbs

  • Test basé sur une distribution normale

  • Test itératif : aussi longtemps que des valeurs aberrantes sont détectées Hypothèses du test:

      Ho : Il n’y a pas d’outlier dans l’échantillon

  H1 : Il y a au moins un outlier dans l’échantillon

Statistique du test : \[G=max_i|X_i-\bar{X}|/\sigma\]

  • Le test de grubbs est disponible dans R : voir la fonction grubbs.test dans le package outliers

  • Le test de grubbs est disponible dans STATA : voir la commande grubbs dans le complément grubbs

Détection- Test de Dixon

  • Méthode adapté pour les petits échantillons (n<=30)

Principe:

Notons tout d’abord qu’il peut s’appliquer aussi bien pour une série statistique à une variable \((x_i)\) que pour une série statistique bivariée \((x_i ; y i)\). Dans le premier cas, les valeurs x_i étant rangées dans l’ordre croissant, le test de Dixon va détecter la (ou les) valeur(s) aberrante(s), aux extrémités de la distribution Comparer la distance entre la valeur suspectée aberrante et une valeur des plus proches, avec la distance entre la valeur suspectée aberrante et une des valeurs les plus éloignées de l’échantillon. Plus le rapport est élevée, plus la valeur suspectée est aberrante

Valeur critique

On se fixe un seuil de risque \(\alpha\). La valeur critique est notée \(r_{1−\alpha}\), elle est définie par : \(P(R \leq r_{1−\alpha}) = 1 − \alpha\) et elle est donnée par la table en fin d’article.

Exemple d’utilisation de la table : \(n = 8\) et \(\alpha = 0,01\). Dans le cas de la recherche d’une valeur aberrante, la table de Dixon indique que pour \(n = 8\) et \(\alpha = 0,01\) , la valeur critique est \(r_{0,99} = 0,59\). Cela signifie que si l’on prélève aléatoirement un échantillon de taille 8 dans une population dans laquelle les données sont distribuées normalement alors la probabilité que \(R\) prenne une valeur inférieure ou égal à 0,59 est 0,99. #### Règle de décision

  • Si \(R_{obs} > r_{1−\alpha}\), on rejette \(H_0\) , donc la valeur suspectée est aberrante.

  • Si \(R_{obs} \leq r_{1−\alpha}\), on n’est pas en mesure de rejeter \(H_0\) .

Statistique du test

la valeur suspectée aberrante est \(x_1\) la valeur suspectée aberrante est \(x_n\)
\[n<=10\] \[R_obs=\frac{x_2-x_1}{x_n-x_1}\] \[R_obs=\frac{x_n-x_{n-1}}{x_n-x_1}\]
\[n>10\] \[R_obs=\frac{x_3-x_1}{x_{n-2}-x_1}\] \[R_obs=\frac{x_n-x_{n-2}}{x_n-x_3}\]

  • Le test de Dixon peut s’intéresser au cas bi-varié.

  • Le test est basé la distribution des résidus de la régression simple de y sur x.

Exemple1

Dans la fabrication de comprimés effervescents, il est prévu que chaque comprimé doit contenir 1 625 mg de bicarbonate de sodium. Afin de contrôler la fabrication de ces médicaments, on a prélevé un échantillon de 10 comprimés et on a mesuré la quantité de bicarbonate de sodium en mg pour chacun d’eux. Les résultats obtenus sont résumés dans le tableau suivant:

1620 1621 1623 1628 1633 1635 1637 1641 1643 1659

On peut demander aux étudiants de réaliser un graphique sur un axe gradué pour détecter quelle(s) valeur(s) semble(nt) aberrante(s).

On effectue un test de Dixon au seuil de risque 0,05 pour tester si la valeur supérieure 1659 est aberrante. On teste les deux hypothèses :

H0 : “1659 n’est pas une valeur aberrante.”

H1 : “1659 est une valeur aberrante.”

\(n = 10\) donc on utilise la variable aléatoire \(R\) qui prend comme valeur observée

\[R_{obs} = \frac{x_n − x_{n − 1}}{x_n − x_1}, soit Robs = R_{obs} = \frac{x_10 − x_9}{x10 − x1}\] qui est égale à 0,410.

D’après la table, la valeur critique est r_{0,95}= 0,412. Comme \(0,41 < 0,412\) : on n’est pas en mesure de rejeter \(H_0\) . La valeur \(1659\) ne peut pas être considérée comme aberrante, au seuil de \(0,05\).

Exemple 2

Une entreprise étudie la possibilité de lancer sur le marché un yaourt à la rhubarbe. Elle réalise des mesures de pH sur un échantillon de 11 pots. Les mesures observées sont les suivantes :

5,40 5,70 6,15 6,16 6,18 6,25 6,43 6,45 6,45 6,60 6,75

Existe-il une valeur aberrante ? Dans un premier temps, nous allons effectuer un test de Dixon au seuil de risque 0,05 sur la valeur \(x_1 = 5,40\) de manière ensuite à justifier la distinction qui doit être faite entre \(n \leq 10\) et \(n > 10\) pour la valeur observée de R.

Le nombre d’observations est ici 11 qui est supérieur à 10, que se passerait-il si on utilisait la valeur observée du cas \(n \leq 10\) ? \(\frac{x_2 − x_1}{x_{11} − x_1}\sim 0,222 (à 10^{-3} près)\)

Bien que nous ne disposions pas de la valeur tabulée pour \(n = 11\), il semble évident que la valeur critique \(r_{0,95}\) serait largement supérieure à 0,222. Il faudrait donc en conclure que 5,40 n’est pas une valeur aberrante.

Cependant, si on élimine cette valeur de l’échantillon et que l’on effectue un test de Dixon au seuil 0,05 sur la valeur \(x_2 = 5,70\) en considérant les 10 valeurs restantes, on observe alors que 5,70 est une valeur aberrante (\(R_{obs} \sim 0,429\) et \(r_0,95 = 0,412\)).

Cette situation invite les étudiants à s’interroger sur cette anomalie car il parait évident que si la deuxième valeur est aberrante, la première l’est tout autant. L’erreur de décision qui est faite en utilisant \(\frac{x_2 − x_1}{x_{11} − x_1}\) se justifie par le fait que les deux premières valeurs sont proches et toutes deux aberrantes. On vérifie alors que l’utilisation de \(\frac{x_3 − x_1}{x_{10} − x_1}\) permet de conclure à l’aberration de la première valeur (\(Robs \sim 0,714\) et \(r_{0,95} = 0,637\)). Pour des échantillons de taille strictement supérieure à 10, le calcul de \(Robs = \frac{x3 − x1}{x_{n − 2} − x_{1}}\) prend en compte la possibilité d’avoir deux valeurs aberrantes inférieures (\(x_1\) et \(x_2\)). Cette situation est plus rare avec des tailles d’échantillon faibles (\(n \leq 10\)).

n/p 0,01 0,05
3 0,988 0,941
4 0,889 0,765
5 0,780 0,642
6 0,698 0,560
7 0,637 0,507
8 0,590 0,468
9 0,555 0,437
10 0,527 0,412
11 0,745 0,637
12 0,704 0,600
13 0,670 0,570
14 0,641 0,546
15 0,616 0,525
16 0,595 0,507
17 0,577 0,490
18 0,561 0,475
19 0,547 0,462
20 0,535 0,450
21 0,524 0,440
22 0,514 0,430
23 0,505 0,421
24 0,497 0,413
25 0,489 0,406
26 0,486 0,399
27 0,475 0,393
28 0,469 0,387
29 0,463 0,381
30 0,457 0,376

Détection- Hampel

considérées comme outliers, les valeurs en dehors de l’intervalle constitué par la médiane, plus ou moins k déviation absolue de médiane

\[I=[mediane−𝑘\times e_𝑚;mediane+𝑘\times𝑒_𝑚]\]

  • \(e_𝑚\): Ecart moyen médian
  • \(k\) est généralement 3

Cas pratiques

library(ggplot2)
## Warning: le package 'ggplot2' a été compilé avec la version R 4.2.3
library(VIM)
library(mice)
library(dplyr)

Import data

data("mpg")

Structure de la base

str(mpg)
## tibble [234 × 11] (S3: tbl_df/tbl/data.frame)
##  $ manufacturer: chr [1:234] "audi" "audi" "audi" "audi" ...
##  $ model       : chr [1:234] "a4" "a4" "a4" "a4" ...
##  $ displ       : num [1:234] 1.8 1.8 2 2 2.8 2.8 3.1 1.8 1.8 2 ...
##  $ year        : int [1:234] 1999 1999 2008 2008 1999 1999 2008 1999 1999 2008 ...
##  $ cyl         : int [1:234] 4 4 4 4 6 6 6 4 4 4 ...
##  $ trans       : chr [1:234] "auto(l5)" "manual(m5)" "manual(m6)" "auto(av)" ...
##  $ drv         : chr [1:234] "f" "f" "f" "f" ...
##  $ cty         : int [1:234] 18 21 20 21 16 18 18 18 16 20 ...
##  $ hwy         : int [1:234] 29 29 31 30 26 26 27 26 25 28 ...
##  $ fl          : chr [1:234] "p" "p" "p" "p" ...
##  $ class       : chr [1:234] "compact" "compact" "compact" "compact" ...
head(force(mpg),10)
## # A tibble: 10 × 11
##    manufacturer model      displ  year   cyl trans drv     cty   hwy fl    class
##    <chr>        <chr>      <dbl> <int> <int> <chr> <chr> <int> <int> <chr> <chr>
##  1 audi         a4           1.8  1999     4 auto… f        18    29 p     comp…
##  2 audi         a4           1.8  1999     4 manu… f        21    29 p     comp…
##  3 audi         a4           2    2008     4 manu… f        20    31 p     comp…
##  4 audi         a4           2    2008     4 auto… f        21    30 p     comp…
##  5 audi         a4           2.8  1999     6 auto… f        16    26 p     comp…
##  6 audi         a4           2.8  1999     6 manu… f        18    26 p     comp…
##  7 audi         a4           3.1  2008     6 auto… f        18    27 p     comp…
##  8 audi         a4 quattro   1.8  1999     4 manu… 4        18    26 p     comp…
##  9 audi         a4 quattro   1.8  1999     4 auto… 4        16    25 p     comp…
## 10 audi         a4 quattro   2    2008     4 manu… 4        20    28 p     comp…

Type de voiture

table(mpg$class)#,mpg$drv)
## 
##    2seater    compact    midsize    minivan     pickup subcompact        suv 
##          5         47         41         11         33         35         62

Méthode graphique

boxplot(mpg$hwy)

Boxplot désagrégé

boxplot(mpg$hwy~mpg$class)

ggplot2 est dans tydiverse

ggplot(mpg,aes(x=class,y=hwy))+
  geom_boxplot()+
  geom_jitter()+
  theme_dark()

Récupération des valeurs atypiques

b <- boxplot(mpg$hwy~mpg$class)

str(b)
## List of 6
##  $ stats: num [1:5, 1:7] 23 24 25 26 26 23 26 27 29 33 ...
##  $ n    : num [1:7] 5 47 41 11 33 35 62
##  $ conf : num [1:2, 1:7] 23.6 26.4 26.3 27.7 26.3 ...
##  $ out  : num [1:19] 35 37 35 44 17 12 12 12 22 44 ...
##  $ group: num [1:19] 2 2 2 2 4 5 5 5 5 6 ...
##  $ names: chr [1:7] "2seater" "compact" "midsize" "minivan" ...
b$stats
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]   23   23   23   21   15 20.0 14.0
## [2,]   24   26   26   22   16 24.5 17.0
## [3,]   25   27   27   23   17 26.0 17.5
## [4,]   26   29   29   24   18 30.5 19.0
## [5,]   26   33   32   24   20 36.0 22.0
b$out
##  [1] 35 37 35 44 17 12 12 12 22 44 41 12 12 25 24 27 25 26 23
head(mpg)
## # A tibble: 6 × 11
##   manufacturer model displ  year   cyl trans      drv     cty   hwy fl    class 
##   <chr>        <chr> <dbl> <int> <int> <chr>      <chr> <int> <int> <chr> <chr> 
## 1 audi         a4      1.8  1999     4 auto(l5)   f        18    29 p     compa…
## 2 audi         a4      1.8  1999     4 manual(m5) f        21    29 p     compa…
## 3 audi         a4      2    2008     4 manual(m6) f        20    31 p     compa…
## 4 audi         a4      2    2008     4 auto(av)   f        21    30 p     compa…
## 5 audi         a4      2.8  1999     6 auto(l5)   f        16    26 p     compa…
## 6 audi         a4      2.8  1999     6 manual(m5) f        18    26 p     compa…

Enlever les atypiques

Methode 1 : les out and group ou Methode 2 : les stats

c <- t(b$stats[c(1,5),])
str(c)
##  num [1:7, 1:2] 23 23 23 21 15 20 14 26 33 32 ...
c <- as.data.frame(c)
c$class <- b$names
c
##   V1 V2      class
## 1 23 26    2seater
## 2 23 33    compact
## 3 23 32    midsize
## 4 21 24    minivan
## 5 15 20     pickup
## 6 20 36 subcompact
## 7 14 22        suv
library("dplyr")
mpg <- mpg%>%left_join(c)
## Joining, by = "class"
mpg <- mpg%>%mutate(outlier=(hwy>V2 | hwy<V1))

Autrement : on construit nous même \(mean\pm 1.5(q3-q1)\)

mpg <- mpg %>% group_by(class) %>% 
  mutate(outlier2=(hwy>quantile(hwy,0.75,na.rm = T)+1.5*IQR(hwy,na.rm = T) | 
                   hwy<quantile(hwy,0.25,na.rm = T)-1.5*IQR(hwy,na.rm = T)))

table(mpg$outlier,mpg$outlier)
##        
##         FALSE TRUE
##   FALSE   215    0
##   TRUE      0   19

Méthode des quantiles

mpg <- mpg %>% group_by(class) %>% 
  mutate(outlier3=(hwy>quantile(hwy,0.975,na.rm = T) | hwy<quantile(hwy,0.025,na.rm = T)))

mpg$outlier3==mpg$outlier
##   [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [13]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [25] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [37]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [49]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [61]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE
##  [73]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [85]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [97]  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [109]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [121]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [133]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [145] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [157]  TRUE  TRUE  TRUE FALSE FALSE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE
## [169]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [181]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [193]  TRUE  TRUE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [205]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [217]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
## [229]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

##Méthode Hampel

mpg <- mpg %>% group_by(class) %>% 
  mutate(outlierH=(hwy>median(hwy,na.rm = T)+3*mad(hwy) | hwy<median(hwy,na.rm = T)-3*mad(hwy)))

Test de GRUBBS : seuls les groups à effectif >30 sont éligibles

Extraction des données de la class

suv <- mpg %>% filter(class=="suv")
nrow(suv)
## [1] 62

Sans itération d’abord

library(outliers)
grubbs.test(suv$hwy)
## 
##  Grubbs test for one outlier
## 
## data:  suv$hwy
## G = 2.97886, U = 0.85215, p-value = 0.06296
## alternative hypothesis: highest value 27 is an outlier
grubbs.test(suv$hwy,opposite = T)
## 
##  Grubbs test for one outlier
## 
## data:  suv$hwy
## G = 2.05812, U = 0.92942, p-value = 1
## alternative hypothesis: lowest value 12 is an outlier

Paramétre type de GRUBBS

  • 10 : une valeur d’un coté
  • 11 : min et max
  • 20 : Deux valeurs
  • opposite = True : on inverse
compact <- mpg %>% filter(class=="compact")
grubbs.test(compact[1:30,]$hwy,type=20)
## 
##  Grubbs test for two outliers
## 
## data:  compact[1:30, ]$hwy
## U = 0.68351, p-value = 0.2143
## alternative hypothesis: highest values 31 , 33 are outliers

Avec itérartion

while (grubbs.test(compact$hwy)$p.value<0.05){
  compact <- compact %>% filter(hwy!=max(hwy))
}
nrow(compact)
## [1] 45
last_tes <- max(compact$hwy)
last_tes
## [1] 35

Ma méthode

while (grubbs.test(compact$hwy)$p.value<0.05){
  compact <- compact %>% filter(hwy!=max(hwy))
  if (grubbs.test(compact$hwy)$p.value>=0.05){
    print(max(compact$hwy))
  }
}

TEST DE DIXON

table(mpg$class)
## 
##    2seater    compact    midsize    minivan     pickup subcompact        suv 
##          5         47         41         11         33         35         62
minivan <- mpg %>% filter(class=="minivan")
dixon.test(minivan$hwy)
## 
##  Dixon test for outliers
## 
## data:  minivan$hwy
## Q = 0.71429, p-value = 0.009622
## alternative hypothesis: lowest value 17 is an outlier
while (dixon.test(minivan$hwy)$p.value<0.05){
  minivan <- minivan %>% filter(hwy!=min(hwy))
}
nrow(minivan)
## [1] 10
last_tes <- max(minivan$hwy)
last_tes
## [1] 24
table(minivan$hwy)
## 
## 21 22 23 24 
##  1  3  2  4