Laboratorio 03: Detección de valores atípicos

Sobre los datos

1.a) Explore y explique qué características posee el dataset MPI_national.csv.

data = read.csv('https://raw.githubusercontent.com/dm-uba/dm-uba.github.io/master/2021/laboratorios/LAB03/MPI_national.csv', sep = ",",header = TRUE)

Observamos la estructura del dataset y las características de las variables numéricas:

library(tidyverse)
library(knitr)

# Structura de variables no numéricas
# Para usar nagate y %>% usar librería tidyverse
data %>% select_if(negate(is.numeric)) %>% str(data)

## 'data.frame':    102 obs. of  2 variables:
##  $ ISO    : chr  "KAZ" "SRB" "KGZ" "TUN" ...
##  $ Country: chr  "Kazakhstan" "Serbia" "Kyrgyzstan" "Tunisia" ...

# Veamos el rango de las variables numéricas
# Kable del paquete knitr sirve para mostrar en formato tabla
# El rango de las variables para zonas rurales es mayor que el rango de la misma variable para zona urbanas. 

kable(data %>% select_if(is.numeric) %>% summarise_each(funs(range)))

MPI.Urban	Headcount.Ratio.Urban	Intensity.of.Deprivation.Urban	MPI.Rural	Headcount.Ratio.Rural	Intensity.of.Deprivation.Rural
0.000	0.0	33.3	0.000	0.09	33.3
0.459	82.5	55.7	0.669	96.92	69.5

# Cantidad de paises
length(unique(data$Country))

## [1] 102

Se cuenta con tres tipos de mediciones recolectadas para zonas urbanas y rurales en 102 países en desarrollo. Existe un solo registro por país (largo de dataset coincide con cantidad de paises). Se desconoce la fecha del relevamiento de los datos. El rango de las variables para zonas rurales es mayor que el rango de la misma variable para zona urbanas.

1.b) Cuál es la distribución de las variables:

# Calculamos el desvío standard y comparamos las zonas para cada medida
sd <- round(apply(data[,c(3:8)],2, sd),2)
sd[order(names(sd))]

##          Headcount.Ratio.Rural          Headcount.Ratio.Urban 
##                          33.27                          18.50 
## Intensity.of.Deprivation.Rural Intensity.of.Deprivation.Urban 
##                           8.78                           5.14 
##                      MPI.Rural                      MPI.Urban 
##                           0.20                           0.09

Las zonas rurales presentan mayor dispersión comparadas con la misma métrica en zonas urbanas.

# Para el histograma usaremos ggplot2
library(ggplot2)

# Instalamos una librería que permite visualizar varios charts en una ventana
# https://github.com/thomasp85/patchwork

# devtools::install_github("thomasp85/patchwork")
library(patchwork)

## Warning: package 'patchwork' was built under R version 4.0.5

# Generamos los histogramas para zonas urbanas
tamaño_texto = 7 
cant_bins = 30

u1 <- qplot(data$MPI.Urban, geom="histogram", main = "Hist MPI.Urban",  
            xlab = "MPI.Urban", ylab = "Frecuencia", binwidth=diff(range(data$MPI.Urban))/cant_bins) + theme(text = element_text(size = tamaño_texto))

u2 <- qplot(data$Headcount.Ratio.Urban, main = "Hist HC Ratio Urban", geom="histogram", xlab = "Headcount Ratio Urban", ylab = "Frecuencia", binwidth=diff(range(data$Headcount.Ratio.Urban))/cant_bins) + theme(text = element_text(size = tamaño_texto))

u3 <- qplot(data$Intensity.of.Deprivation.Urban, main = "Hist. Intensity of Depr. Urban", geom="histogram",  xlab = "Intensity of Depr. Urban", ylab = "Frecuencia",binwidth=diff(range(data$Intensity.of.Deprivation.Urban))/cant_bins) + theme(text = element_text(size = tamaño_texto))

# Y los mismos para zonas Rurales
r1 <- qplot(data$MPI.Rural, geom="histogram", main = "Hist MPI Rural",  
            xlab = "MPI.Rural", ylab = "Frecuencia", binwidth=diff(range(data$MPI.Rural))/cant_bins) + theme(text = element_text(size = tamaño_texto))

r2 <- qplot(data$Headcount.Ratio.Rural, main = "Hist HC Ratio Rural", geom="histogram", xlab = "Headcount Ratio Rural", ylab = "Frecuencia", binwidth=diff(range(data$Headcount.Ratio.Rural))/cant_bins) + theme(text = element_text(size = tamaño_texto))

r3 <- qplot(data$Intensity.of.Deprivation.Rural, main = "Hist. Int. of Depr. Rural", geom="histogram",  xlab = "Intensity of Depr. Rural", ylab = "Frecuencia",binwidth=diff(range(data$Intensity.of.Deprivation.Rural))/cant_bins) + theme(text = element_text(size = tamaño_texto))


# Y los graficamos en paralelo para facilitar la comparación
(u1 | u2 | u3) /
  (r1 | r2 | r3)

Las zonas urbanas presentan una asimetría a derecha más marcada que las zonas rurales para MPI y Ratio de pobreza per capita, mientras que la distribución de Intensidad de Privaciones para ambas zonas se asemeja.

Tratamiento de outliers

2.a) ¿Existen atributos que poseen valores atípicos?

Normalizamos de variables con z-score para que sean comparables. Este método también nos informa sobre a cuántos desvíos de la media se distribuyen los datos.

for(i in 3:ncol(data)) {
data[,i] <- (data[,i]-mean(data[,i]))/(sd(data[,i]))}

# Configuramos el tamaño del lienzo (ver "par" del paquete graphics)
par(mfrow=c(1, 1), mar=c(5, 12, 4, 2) ) 

# Y graficamos un solo boxplot con todas las medidas
boxplot(data[,c(3:8)],horizontal = TRUE, las=1, xlab="z-score")
stripchart(data[,c(3:8)], vertical = FALSE ,
    method = "jitter", add = TRUE, pch = 20, col = 'blue')

# Que paises son outliers segun boxplot de MPI Urban
MPI.Urban_Umbral = boxplot(data[,c(3)])$stats[5]

data["Country"][data$MPI.Urban > MPI.Urban_Umbral,]

## [1] "Chad"        "South Sudan"

# Que paises son outliers segun boxplot de Headcount Ratio Urban
HC.Ratio.Urban_umbral = boxplot(data[,c(4)])$stats[5]

data["Country"][data$Headcount.Ratio.Urban > HC.Ratio.Urban_umbral,]

## [1] "South Sudan"

Según el boxplot para MPI Urban Chad y South Sudan son países atípicos, mientras que si observamos el ratio de pobreza per capita para zonas urbanas (HC Ratio Urban), el país considerado outlier es South Sudan.

2.b) Seleccione una variable con outliers y analice utilizando los métodos univariados IRQ, SD, y Z-score

### Focalicemos el análisis en MPI Urban ###
variable = data$MPI.Urban

# Rango intercuartil
Q1 = as.numeric(quantile(variable)["25%"])
Q3 = as.numeric(quantile(variable))[4]

IQR = Q3 - Q1

lim_sup = Q3 + 1.5*IQR
lim_inf = Q1 - 1.5*IQR

print(lim_sup)

## [1] 2.403141

# Veamos de qué paises estamos hablando
# Coinciden con los detectados en el boxplot porque ambos usan 1.5 * IQR

sort(data$Country[variable > lim_sup]) 

## [1] "Chad"        "South Sudan"

# con Z Score
# Ya hicimos el escalado a z-score

umbral_zscore =3
unique(sort(variable[variable> umbral_zscore])) 

## [1] 4.062822

# Veamos de qué país estamos hablando
sort(data$Country[variable > umbral_zscore]) 

## [1] "South Sudan"

Utilizando z-score el pais que representa un valor atipico es South Sudan

# Bonus Track
# veamos que pasa con z-score modificado

data_original = read.csv('https://raw.githubusercontent.com/dm-uba/dm-uba.github.io/master/2021/laboratorios/LAB03/MPI_national.csv', sep = ",",header = TRUE)

variable_zm = data_original$MPI.Urban

mediana = median(variable_zm)
MAD = median(abs(variable_zm-mediana))
Zm = ((0.6745*(variable_zm - median(variable_zm))) / MAD)

umbral_Zm = 3.5

# Cuantos outliers identifica?
length(Zm[Zm > umbral_Zm])

## [1] 11

# Que paises conforman la lista?
sort(data$Country[Zm > umbral_Zm]) 

##  [1] "Burkina Faso"                      "Burundi"                          
##  [3] "Central African Republic"          "Chad"                             
##  [5] "Congo, Democratic Republic of the" "Ethiopia"                         
##  [7] "Liberia"                           "Niger"                            
##  [9] "Sierra Leone"                      "Somalia"                          
## [11] "South Sudan"

# Observemos el diagrama de dispersión

plot(Zm, main="Diag. de Dispersión Z Score Modificado")
abline(h=c(mediana, umbral_Zm), col=c("red","green"),lty=5)
legend(1,8 , legend=c("Mediana", "Zm Umbral"),col=c("red","green"), lty=5)

2.c) Observe qué ocurre con la distribución de la feature elegida en caso de eliminar los outliers. Grafique un boxplot de con la nueva distribución. Concluya al respecto

# Observemos la distribución de MPI Urban quitando los outliers identificados en el boxplot
MPI.Urban_nuevo_bp = boxplot(data[c(3)][data[c(3)] <= MPI.Urban_Umbral],main="MPI Urban sin outliers", ylab="z-score", xlab="MPI Urban acotada") 

# Luego de eliminar los outliers, el nuevo valor atipico es:
MPI.Urban_nuevo_bp$out

## [1] 2.291073

# Este nuevo outlier corresponde a:
data["Country"][data[, "MPI.Urban"] == MPI.Urban_nuevo_bp$out,]

## [1] "Somalia"

Dada la asimetría de la distribución, luego de eliminar outliers, Somalía es ahora clasificado como atípico.

2.d) Extienda el análisis a 3 variables y analice si existen valores atípicos utilizando algún método multivariado.

# Distancia Mahalanobis para zonas urbanas
vector_medias = colMeans(data[,c(3:5)]) 
matriz_var_cov = cov(data[,c(3:5)])

# Creamos una variable con la distancia
data$maha = sqrt(mahalanobis(data[,c(3:5)],vector_medias,matriz_var_cov))

# Los 3 registros mas distantes
top_maha <- head(data[order(data$maha,decreasing = TRUE),],3)
print(top_maha$Country)

## [1] "South Sudan" "Chad"        "Montenegro"

¿Qué países son considerados outliers según LOF?

library(Rlof)

data$LOF_score<-lof(data[,c(3:5)], k=5)
top_LOF <- head(data[order(data$LOF_score,decreasing = TRUE),],3)
print(top_LOF$Country)

## [1] "South Sudan"         "Montenegro"          "Trinidad and Tobago"

Veamos donde se ubican los puntos

library(scatterplot3d)

lof_3d <- scatterplot3d(x = data[,3], y = data[,4], z = data[,5],
color=ifelse(data$LOF_score>=min(top_LOF$LOF_score),"red","black"), xlab="z-score MPI Urban",ylab="z-score HC Ratio Urban",zlab="z-score Intensity of Deprivation Urban" , main="3D Scatterplot LOF Scores")

# Uso el método xyz.convert para convertir el objeto a 2 variables
lof_3d.coords <- lof_3d$xyz.convert(data[,3], data[,4], data[,5]) 

# Asigno como etiquetas los nombres de los países
text(lof_3d.coords$x, 
     lof_3d.coords$y,             
     labels = ifelse(data$LOF_score>=min(top_LOF$LOF_score),data$Country,""),              
     cex = .5, 
     pos = 4)  

Ahora hagamos el análisis utilizando Isolation Forest

library(isotree)
library(ggplot2)

# Entrenamos el modelo IForest y nos quedamos con el score
# Utilizamos el conjunto de datos en su totalidad

data$iforest= isolation.forest(data[,c(3:5)], ntrees = 100, ndim=1, random_seed = 13, output_score=TRUE)$score

library(dplyr)
library(highcharter)

# Set highcharter options
options(highcharter.theme = hc_theme_smpl(tooltip = list(valueDecimals = 2)))

hchart(data$iforest, name="iForest score", color="green") %>%
  hc_title(text = "Histograma Scores Isolation Forest")

top_iforest <- head(data[order(data$iforest,decreasing = TRUE),],3)
print(top_iforest$Country)

## [1] "South Sudan" "Chad"        "Somalia"

# 
iforest_3d <- scatterplot3d(x = data[,3], y = data[,4], z = data[,5],
color=ifelse(data$iforest>=min(top_iforest$iforest),"red","black"), xlab="z-score MPI Urban",ylab="z-score HC Ratio Urban",zlab="z-score Intensity of Deprivation Urban", main="3D Scatterplot Isolation Forest Scores")

# Uso el método xyz.convert para convertir el objeto a 2 variables
iforest_3d.coords <- iforest_3d$xyz.convert(data[,3], data[,4], data[,5]) 

# Asigno el nombre de país como etiquetas
text(iforest_3d.coords$x, 
     iforest_3d.coords$y,             
     labels = ifelse(data$iforest>=min(top_iforest$iforest),data$Country,""),              
     cex = .5, 
     pos = 4) 

## Warning: package 'rgl' was built under R version 4.0.5

# Ejemplo de un gráfico 3D interactivo
# Prueben el código como R Script 

library(rgl)

data$Outlier_iforest = ifelse(data$iforest>=min(top_iforest$iforest),"S","N" )

colors <- c("royalblue1", "red")
data$color <- colors[ as.numeric( as.factor(data$Outlier_iforest) ) ]

# Para poder interactuar uso plot3d
plot3d(x=data[,3], y=data[,4],z=data[,5], col = data$color, type = "s", size = 1, xlab="z-score MPI Urban",ylab="z-score HC Ratio Urban",zlab="z-score Intensity of Deprivation Urban" )

# rglwidget()

Veamos ahora cómo sería la implementación utilizando un muestreo en las iteraciones de los árboles. No podemos utilizar el tamaño de muestra recomendado (256) porque solo contamos con 102 registros en el conjunto de datos:

#Ajusto el modelo
# Agregamos el parámetro sample_size

iforest_sample = isolation.forest(data[,c(3:5)], ntrees = 100, ndim=1, random_seed = 13, sample_size = 50)

# Utilizo el método predict para calcular el score sobre todos los registros
data$iforest_pred = predict(iforest_sample, data[,c(3:5)])

# Qué países identifico?
top_iforest_sample <- head(data[order(data$iforest_pred,decreasing = TRUE),],3)
print(top_iforest_sample$Country)

## [1] "South Sudan" "Chad"        "Somalia"