R Markdown

La dedection d’outliers est une branche particulière dans la statistique, bien que souvent négligé lors des analyses de données, les méthodes d’outliers servent à identifier les variables qui ne correspondent pas à la majorité des données mais qui cependant peuvent avoir un effet important sur certaines analyses.

Durant cette introduction nous commencerons par étudier les principes de bases de la détection d’outliers. Nous ferons des travaux pratiques concernant la détection d’outliers dans le cas d’une variable et à 2 variables quantitatives.

Définition : Une observation est dîtes outlier si elle dévie de manière significative des autres observations. On suspecte donc que cette observation a été générée par un autre mécanisme que les autres observations.

Remarque : Un outlier n’a pas forcément d’influence sur les données tout comme une observation avec une grande influence n’est pas forcément un outlier.

L’étude des outliers a deux objectifs ; Un premier qui consiste à définir des méthodes statistiques qui ne soient pas trop sensibles aux valeurs aberrantes. Le second objectif est de détecter les valeurs aberrantes afin de les analyser comme il peut être le cas dans la détection des maladies, de la fraude ou des produits défectueux.

L’objectif sera d’apprendre les principes théoriques en même temps que la partie pratique.

#Lobjectif de cette première introduction est d'étudier les outliers sur une variable. ici les données sont le nombre de publications de différents instituts.

setwd("C:/Users/Vicelord/Downloads")

#------------------------------------------------------------------------------------------------------------#
#                                           Univariate approach                                              #
#------------------------------------------------------------------------------------------------------------#
#------------------------------------------------------------------------------------------------------------#
#                                         1.3 Description de la base de données                                 #
#------------------------------------------------------------------------------------------------------------#
# Import the data
data <- read.csv("Scimago_data.csv",row.names=1)

# Description des données
summary(data)
##   Publication           IC              Q1              NI       
##  Min.   : 26035   Min.   :10.40   Min.   :17.90   Min.   :0.500  
##  1st Qu.: 29626   1st Qu.:24.38   1st Qu.:55.42   1st Qu.:1.200  
##  Median : 33756   Median :27.30   Median :67.95   Median :1.800  
##  Mean   : 40204   Mean   :30.43   Mean   :61.00   Mean   :1.618  
##  3rd Qu.: 40182   3rd Qu.:34.90   3rd Qu.:71.95   3rd Qu.:2.000  
##  Max.   :144269   Max.   :65.00   Max.   :84.30   Max.   :2.600  
##       Exc       
##  Min.   : 3.10  
##  1st Qu.:17.75  
##  Median :24.15  
##  Mean   :22.09  
##  3rd Qu.:28.45  
##  Max.   :40.10
## on remarque que le Max de publication est largement supérieur à la médiane et au 3 eme quartile.

boxplot(data, main = "Description of the data (Scimago 2011)")

# Données standardisées
data.std = scale(data, center = TRUE, scale = TRUE)
summary(data.std)
##   Publication              IC                Q1                NI         
##  Min.   :-0.617998   Min.   :-1.7792   Min.   :-2.5573   Min.   :-2.0157  
##  1st Qu.:-0.461386   1st Qu.:-0.5381   1st Qu.:-0.3308   1st Qu.:-0.7536  
##  Median :-0.281244   Median :-0.2783   Median : 0.4124   Median : 0.3281  
##  Mean   : 0.000000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.:-0.000973   3rd Qu.: 0.3966   3rd Qu.: 0.6497   3rd Qu.: 0.6887  
##  Max.   : 4.538808   Max.   : 3.0698   Max.   : 1.3825   Max.   : 1.7705  
##       Exc         
##  Min.   :-2.1215  
##  1st Qu.:-0.4849  
##  Median : 0.2301  
##  Mean   : 0.0000  
##  3rd Qu.: 0.7105  
##  Max.   : 2.0121
boxplot(data.std)

head(data)
##                             Publication   IC   Q1  NI  Exc
## Chinese Academy of Sciences      144269 21.5 40.5 0.9 11.3
## CNRS                             130977 49.0 61.9 1.4 18.7
## Russian Academy of Sciences       88907 35.0 24.2 0.5  5.9
## Harvard University                69995 34.4 79.0 2.4 35.7
## Inst_5                            49987 65.0 72.2 1.8 29.3
## Inst_6                            48947 26.3 56.7 1.2 17.9
# Les observations sont ordonnées, on supprime donc les 4 premieres observations.
dataclean = data[-c(1:4),]
boxplot(scale(dataclean, center = TRUE, scale = TRUE))

# Nous avons plus d'outliers sur notre première variable publication.

Location estimator et Courbe de sensibilité

Une notion très importante en outlier est la courbe de sensibilité, elle permet notamment d’identifier l’évolution de certaines valeurs statistiques comme la moyenne par la présence de nouvelle observation. Dans cette partie nous allons introduire le principe de “location estimator”.

Presnetons les 3 méthodes de location estimators :

Prenons un exemple simple pour comprendre la différence entre la moyenne, le alpha trimmed et la Winsorized mean Considérons l’échantillon : 2, 4, 5, 10, 200.

Then for ?? = 0.2: mean = 36.83 Comme on a fixé alpha à 0.2, pour la trimmed mean on conserve seulement 80% de l’échantillon en enlevant les 10% les plus faible et les 10% les plus fortes. trimmed mean = (4+5+10)/3 = 6.33 Pour la Winsorized mean le principe est différent. Elle consiste à remplacer les 10% les plus faible (forte) par la valeur étant le quantile à 10 % (90%) de l’échantillon. winsorized mean = (4+4+5+10+10)/5 = 6.6

Ces différentes méthodes permettent le calcul de location plus robuste que la moyenne

#------------------------------------------------------------------------------------------------------------#
#                                          1.4 Sensitivity curve                                             #
#------------------------------------------------------------------------------------------------------------#
#----------------------------------#
#   1.4.1 Location estimator       #

#----------------------------------#
# 1. Les 4 principaux "location estimators"" sont :
# La moyenne, la mediane ainsi que l'alpha trimmed, winsorized,

# Analysons les 3 premières méthodes

# 2. Prenons par exemple 3 nouvelles observations. Une très faible, une correspondant à la mediane et une dernière très elevé

new.obs = c(15000, median(dataclean[,1]),55000)

# Pour chacune des nouvelles observations on va calculer la moyenne, la médiane et la trimmed mean.

Tn1 = t(sapply(new.obs, function(x){
  # to create a new 'Publication' vector containing the new point
  y = c(x,dataclean[,1])
  # to compute the three location estimates when we add the point x
  c(Mean = mean(y), Median = median(y), alphaMean = mean(y, trim = 0.1))
} ))

# The values of the location estimates on the initial 'Publication' variable
Tn = c(Mean = mean(dataclean[,1]), Median = median(dataclean[,1]), alphaMean = mean(dataclean[,1], trim = 0.1))
print(Tn)
##      Mean    Median alphaMean 
##  34262.33  32627.00  33660.71
## 


# Afin de comparer la courbe de sensibillité de chacun , on retourne dans une fonction la fonction qui permet 
# le calcul de la courbe de sensitivité est donné par sc.all.
SC.all = vector()
for (i in 1:length(Tn)){
  SC.all = cbind(SC.all, (nrow(dataclean)+1)*(Tn1[,i]-Tn[i]))
}
colnames(SC.all) = colnames(Tn1)
SC.all
##            Mean Median alphaMean
## [1,] -19262.326 -12831 -8358.420
## [2,]  -1635.326      0 -1245.754
## [3,]  20737.674  12831 14594.451
# On plot la courbe de senibillité 
matplot(new.obs, SC.all, type='l',lty=1:3,col=1:3,xlab="Publication",ylab="Sensitivity curve")
legend("topleft", c("mean","median","trimmed mean"), col=1:3, lty=1:3, lwd=2)

## maintenant à vous d'appliquer , appliquer cette méthode en choissisant 3 valeurs de votre choix et regarder comment évolue la courbe de sensibillité.

######
a=15000
b=55000
c=1000

#####

new.obs = sort(c(seq(a, b,c), median(dataclean$Publication)))

# Compute for each of these three points, the value of the three location estimates when we add the new
# point to the initial 'Publication' variable
Tn1 = t(sapply(new.obs, function(x){
  y = c(x,dataclean[,1])
  c(Mean = mean(y), Median = median(y), alphaMean = mean(y, trim = 0.1))
} ))

# The values of the location estimates on the initial 'Publication' variable
Tn = c(Mean = mean(dataclean[,1]), Median = median(dataclean[,1]), alphaMean = mean(dataclean[,1], trim = 0.1))

# To compute the value of the Sensitivity Curve (SC) at each point for the three location estimators
SC.all = vector()
for (i in 1:length(Tn)){
  SC.all = cbind(SC.all, (nrow(dataclean)+1)*(Tn1[,i]-Tn[i]))
}
colnames(SC.all) = colnames(Tn1)

# Plot the three sensitivity curves on the same plot
matplot(new.obs, SC.all, type='l',lty=1:3,col=1:3,xlab="Publication",ylab="Sensitivity curve")
legend("topleft", c("mean","median","trimmed mean"), col=1:3, lty=1:3, lwd=2)

## Quelle est votre remarque ?


# 5. Quels estimateurs sont robuste?
# Rappelons la definition : un estimateur est robust si et seulement si il est "bounded" c'est à dire délimité par de sbornes fermé
#d'après les différents plot, on remarque que la moyenne n'est pas robuste et est donc sensible aux outliers contrairement aux deux autre estimateurs


# Remarque 
# Que signifie avoir un estimator à 0 ?
# ==> que si nous remplaçons une observation par la valeur correspondante des publications, alors l'estimation de la localisation est la même
# wPour quelles valeurs c'est le cas ?
matplot(new.obs, SC.all, type='l',lty=1:3,col=1:3,xlab="Publication",ylab="Sensitivity curve")
legend("topleft", c("mean","median","trimmed mean"), col=1:3, lty=1:3, lwd=2)
abline( h=0, col="grey", lty=4)
abline(v=Tn, col="grey", lty=4)

Scale estimators

Nous allons maintenant nous intérréser au “scale estimators”. Dans la famille de ces estimateurs, on peut notamment noter l’ecart-type, l’écart interquartile ainsi que la median absolue deviation (MAD) :

MAD = median(|Xi - median(X)|)

#----------------------------------#
# 1.4.2  Scale estimator           #

#----------------------------------#
## Meme exercice que précédamment, on prend 3 nouvelles observations et on regarde comment les mesures d'échelles varient selon leurs courbes de sensibillités

new.obs = c(15000, median(dataclean[,1]),55000)

# On ajoute dans la fonction Sn1 les 3 estimateurs
Sn1 = t(sapply(new.obs, function(x){
  y = c(x,dataclean[,1])
  c(Sd = sd(y), MAD = mad(y), IQR  = IQR(y)/1.349)
} ))
Sn = c(Sd = sd(dataclean[,1]),  MAD = mad(dataclean[,1]), IQR  = IQR(dataclean[,1])/1.349)
SC.all = vector()
for (i in 1:length(Sn)){
  SC.all = cbind(SC.all, (nrow(dataclean)+1)*(Sn1[,i]-Sn[i]))
}
colnames(SC.all) = colnames(Sn1)

matplot(new.obs, SC.all, type='l',lty=1:3,col=1:3,xlab="Publication",ylab="Sensitivity curve")
legend("topleft", c("sd","mad","IQR"), col=1:3, lty=1:3, lwd=2)

# Pareil que précédamment choissisez 3 nouveaux points
# a = ?
#b=?
# c= ?
new.obs = sort(c(seq(a, b, c), median(dataclean$Publication)))

Sn1 = t(sapply(new.obs, function(x){
  y = c(x,dataclean[,1])
  c(Sd = sd(y), MAD = mad(y), IQR  = IQR(y)/1.349)
} ))
Sn = c(Sd = sd(dataclean[,1]),  MAD = mad(dataclean[,1]), IQR  = IQR(dataclean[,1])/1.349)


SC.all = vector()
for (i in 1:length(Sn)){
  SC.all = cbind(SC.all, (nrow(dataclean)+1)*(Sn1[,i]-Sn[i]))
}
colnames(SC.all) = colnames(Sn1)

# Plot the three sensitivity curves on the same plot
matplot(new.obs, SC.all, type='l',lty=1:3,col=1:3,xlab="Publication",ylab="Sensitivity curve")
legend("topleft", c("sd","mad","IQR"), col=1:3, lty=1:3, lwd=2)

# Vos commentaires ?

# 5. Quels estimateurs sont robust ? Même définition que précédamment, on regarde les estimateurs qui sont bounded
# sd: convex function (parabolic)- unbouded 
# mad and IQR: bounded ==> B-robust

Après avoir travaillé sur la detection d’outliers pour une variable ainsi que d’avoir évalué quels étaient les estimateurs robustes en précésence de valeurs abérrantes Nous allons nous intérréssés à un cas plus complexe : le cas détudes multivariées.

Outliers avec données multivariées.

Nous allons appliquer l’essentiel des techniques qui consistes à detceter les outliers sur un plan à deux dimensions.

#------------------------------------------------------------------------------------------------------------#
#                                       2 Multivariate approach                                              #
#------------------------------------------------------------------------------------------------------------#
# Nous allons uttiliser une nouvelle base de donnée appelée ellipse

X = read.csv("data_ellipse.csv")
# On suprime la premiere colonne qui est l'indice des lignes
# Nous avons deux variables V4 et V6
X = X[,-1]
head(X)
##         V4       V6
## 1 26.93385 24.38764
## 2 26.54896 24.94254
## 3 26.30215 23.73536
## 4 26.78286 24.89589
## 5 26.75614 24.33472
## 6 26.51074 24.33692
# Library
library(rrcov)
## Warning: package 'rrcov' was built under R version 3.3.3
## Loading required package: robustbase
## Warning: package 'robustbase' was built under R version 3.3.3
## Scalable Robust Estimators with High Breakdown Point (version 1.4-3)
# Que pouvez vous dire à propos d'outliers ?
# Visiblement les différentes variable ssont bien compacte
summary(X)
##        V4              V6       
##  Min.   :23.20   Min.   :22.30  
##  1st Qu.:25.45   1st Qu.:23.90  
##  Median :26.30   Median :24.70  
##  Mean   :26.01   Mean   :24.67  
##  3rd Qu.:26.77   3rd Qu.:25.40  
##  Max.   :28.00   Max.   :27.40
boxplot(X)

# 2. On plot les scatterplots des 2 variables, que pouvez vous déduire ? 
pairs(X)

plot(X[,1], X[,2], xlim = c(22.5,29.5), ylim = c(21,29))


# Visiblement certaines observations s'éloignent de la masse, mais lesquels peut-on considérer comme outliers ??

#----------------------------------#  
# Graphical representation of  #
# location and scatter estimators  # 
#----------------------------------#
#1. On représnte grapiquement la moyenne des  variables sur un plan à deux dimensions
mu = colMeans(X)
plot(X[,1], X[,2], xlim = c(22.5,29.5), ylim = c(21,29))
points(mu[1], mu[2], pch = 19, col = "red", cex=2)


# Nous allons étudier les outliers à l'aide d'une ellipse

#2. La fonction ci dessous est une fonction générique complexe, elle permet de compacter les observations dans un intervalle et se base selon deux critères :
# la location qui est généralement la moyenne des deux variables sur un plan à deux dimensions, le second paramètre est la covariance c'est à dire la variance #entre les deux variables.
Sigma = var(X)

# GRâce au package rcov on peut avoir :
ellips <- function(loc, cov) {
  dist <- sqrt(qchisq(0.975, 2))
  A <- solve(cov)
  eA <- eigen(A)
  ev <- eA$values
  lambda1 <- max(ev)
  lambda2 <- min(ev)
  eigvect <- eA$vectors[, order(ev)[2]]
  z <- seq(0, 2 * pi, by = 0.01)
  z1 <- dist/sqrt(lambda1) * cos(z)
  z2 <- dist/sqrt(lambda2) * sin(z)
  alfa <- atan(eigvect[2]/eigvect[1])
  r <- matrix(c(cos(alfa), -sin(alfa), sin(alfa), cos(alfa)), 
              ncol = 2)
  t(loc + t(cbind(z1, z2) %*% r))
}

# representation
z1 <- ellips(loc = mu, cov = Sigma)
points(z1, type = "l")
# Interpretation: Sous l'hypothèse gaussienne, on suspecte que 95% des observations sont "non outliers"


#----------------------------------#
# 2.5 The Reweighted Minimum       #
#  Covariance Determinant (RMCD)   #
#----------------------------------#

#4. RMCD method
RMCD = CovMcd(X, alpha = 0.75)
str(RMCD)
## Formal class 'CovMcd' [package "rrcov"] with 22 slots
##   ..@ alpha      : num 0.75
##   ..@ quan       : num 84
##   ..@ best       : int [1:84] 2 4 22 28 29 30 31 32 33 34 ...
##   ..@ raw.cov    : num [1:2, 1:2] 2 1.91 1.91 1.96
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "V4" "V6"
##   .. .. ..$ : chr [1:2] "V4" "V6"
##   ..@ raw.center : Named num [1:2] 25.9 24.7
##   .. ..- attr(*, "names")= chr [1:2] "V4" "V6"
##   ..@ raw.mah    : num [1:111] 14.95 1.75 15.52 4.35 12.42 ...
##   ..@ raw.wt     : num [1:111] 0 1 0 1 0 0 0 0 0 0 ...
##   ..@ raw.cnp2   : num [1:2] 1.83 1.03
##   ..@ cnp2       : num [1:2] 1.73 1.01
##   ..@ iter       : num 500
##   ..@ crit       : num -2.62
##   ..@ wt         : num [1:111] 0 1 0 1 0 1 0 0 0 0 ...
##   ..@ call       : language CovMcd(x = X, alpha = 0.75)
##   ..@ cov        : num [1:2, 1:2] 1.81 1.71 1.71 1.76
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "V4" "V6"
##   .. .. ..$ : chr [1:2] "V4" "V6"
##   ..@ center     : Named num [1:2] 25.9 24.7
##   .. ..- attr(*, "names")= chr [1:2] "V4" "V6"
##   ..@ det        : num -1
##   ..@ n.obs      : int 111
##   ..@ mah        : num [1:111] 12.05 1.33 12.56 3.41 9.97 ...
##   ..@ flag       : NULL
##   ..@ method     : chr "Fast MCD(alpha=0.75 ==> h=84); nsamp = 500; (n,k)mini = (300,5)"
##   ..@ singularity: NULL
##   ..@ X          : num [1:111, 1:2] 26.9 26.5 26.3 26.8 26.8 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:111] "1" "2" "3" "4" ...
##   .. .. ..$ : NULL
z2.RMCD25 <- ellips(loc = RMCD@center, cov = RMCD@cov)
points(z2.RMCD25, type = "l", col="orange")

#5. Test different values of alpha and plot the associated ellipses on the same plot.
plot(X[,1], X[,2], xlim = c(22.5,29.5), ylim = c(21,29), xlab = "X1", ylab = "X2")
points(mu[1], mu[2], pch = 19, col = "red", cex=2)
z1 <- ellips(loc = mu, cov = Sigma)
points(z1, type = "l")


alpha.all = c(0.5, 0.75,0.85,0.95)
col.all = topo.colors(length(alpha.all))
for (i in 1:length(alpha.all) ){
  RMCD = CovMcd(X, alpha = alpha.all[i])
  #print(det(RMCD@cov))
  z <- ellips(loc = RMCD@center, cov = RMCD@cov)
  points(z, type = "l", col = col.all[i])
}
legend("topleft",  legend = c("classic",alpha.all), col = c(1,col.all), lwd=2)