Más sobre ciencia de datos en cienciadedatos.net
Versión PDF
La detección de anomalías (outliers) con PCA y Autoencoders es una estrategia no supervisada para identificar anomalías cuando los datos no están etiquetados, es decir, no se conoce la clasificación real (anomalía - no anomalía) de las observaciones.
Si bien esta estrategia hace uso de PCA o Autoencoders, no utiliza directamente su resultado como forma de detectar anomalías, sino que emplea el error de reconstrucción producido al revertir la reducción de dimensionalidad. El error de reconstrucción como estrategia para detectar anomalías se basa en la siguiente idea: los métodos de reducción de dimensionalidad permiten proyectar las observaciones en un espacio de menor dimensión que el espacio original, a la vez que tratan de conservar la mayor información posible. La forma en que consiguen minimizar la perdida global de información es buscando un nuevo espacio en el que la mayoría de observaciones puedan ser bien representadas.
Los métodos de PCA y Autoencoders crean una función que mapea la posición que ocupa cada observación en el espacio original con el que ocupa en el nuevo espacio generado. Este mapeo funciona en ambas direcciones, por lo que también se puede ir desde el nuevo espacio al espacio original. Solo aquellas observaciones que hayan sido bien proyectadas podrán volver a la posición que ocupaban en el espacio original con una precisión elevada.
Dado que la búsqueda de ese nuevo espacio ha sido guiada por la mayoría de las observaciones, serán las observaciones más próximas al promedio las que mejor puedan ser proyectadas y en consecuencia mejor reconstruidas. Las observaciones anómalas, por el contrario, serán mal proyectadas y su reconstrucción será peor. Es este error de reconstrucción (elevado al cuadrado) el que puede emplearse para identificar anomalías.
Los siguientes packages se emplean a lo largo del documento.
library(R.matlab) # Lectura de archivos .mat
library(rrcov) # PCA robusto
library(factoextra) # Visualización de PCA
library(h2o) # Entrenar autoencoders
library(tidyverse) # Preparación de datos y gráficos
library(progress) # Barra de progresso en loop for
Los datos empleados en este documento se han obtenido de Outlier Detection DataSets (ODDS), un repositorio con sets de datos comúnmente empleados para comparar la capacidad que tienen diferentes algoritmos a la hora de identificar outliers. Shebuti Rayana (2016). ODDS Library [http://odds.cs.stonybrook.edu]. Stony Brook, NY: Stony Brook University, Department of Computer Science.
Todos estos data sets están etiquetados, se conoce si las observaciones son o no anomalías (variable y). Aunque los métodos que se describen en el documento son no supervisados, es decir, no hacen uso de la variable respuesta, conocer la verdadera clasificación permite evaluar su capacidad para identificar correctamente las anomalías.
Los datos están disponibles en formato matlab (.mat). Para leer su contenido se emplea la función readMat()
del paquete R.matlab v3.6.2
.
cardio_mat <- readMat("./datos/cardio.mat")
df_cardio <- as.data.frame(cardio_mat$X)
df_cardio$y <- as.character(cardio_mat$y)
speech_mat <- readMat("./datos/speech.mat")
df_speech <- as.data.frame(speech_mat$X)
df_speech$y <- as.character(speech_mat$y)
shuttle_mat <- readMat("./datos/shuttle.mat")
df_shuttle <- as.data.frame(shuttle_mat$X)
df_shuttle$y <- as.character(shuttle_mat$y)
Principal Component Analysis (PCA) es un método estadístico que permite simplificar la complejidad de espacios muestrales con muchas dimensiones a la vez que conserva su información. Supóngase que existe una muestra con \(n\) individuos cada uno con \(p\) variables (\(X_1\), \(X_2\), …, \(X_p\)), es decir, el espacio muestral tiene \(p\) dimensiones. PCA permite encontrar un número de factores subyacentes, \((z < p)\) que explican aproximadamente lo mismo que las \(p\) variables originales. Donde antes se necesitaban \(p\) variables para caracterizar a cada individuo, ahora bastan \(z\). Estas nuevas variables se llaman componentes principales y tienen la característica de no estar linealmente correlacionadas entre ellas.
Cada componente principal (\(Z_i\)) se obtiene por combinación lineal de las variables originales. Se pueden entender como nuevas variables obtenidas al combinar de una determinada forma las variables originales. La primera componente principal de un grupo de variables (\(X_1\), \(X_2\), …, \(X_p\)) es la combinación lineal normalizada de dichas variables que tiene mayor varianza. Una vez calculada la primera componente (\(Z_1\)), se calcula la segunda (\(Z_2\)) repitiendo el mismo proceso, pero añadiendo la condición de que la combinación lineal no pude estar correlacionada con la primera componente. El proceso se repite de forma iterativa hasta calcular todas las posibles componentes (min(n-1, p)) o hasta que se decida detener el proceso.
La principal limitación que tiene el PCA como método de reducción de dimensionalidad es que solo contempla combinaciones lineales de las variables originales, lo que significa que no es capaz de capturar otro tipo de relaciones.
Para más información sobre PCA consultar Análisis de Componentes Principales (Principal Component Analysis, PCA) y t-SNE
Con la función prcomp()
del paquete stats
se pueden calcular las componentes principales de un dataframe o matriz. El objeto prcomp
almacena la información de los eigenvectors, eigenvalues y la matriz de rotación. Con el método predict()
se pueden proyectar nuevas observaciones.
Tras obtener el PCA, se evalúa su capacidad para reducir la dimensionalidad de los datos explorando los eigenvectors, eigenvalues y la varianza capturada en el proceso.
Información sobre los eigenvalues/varianza
fviz_eig(
X = pca,
choice = "variance",
addlabels = TRUE,
ncp = 21,
ylim = c(0, 100),
main = "Porcentaje de varianza explicada por componente",
xlab = "Componente",
ylab = "Porcentaje de varianza explicada"
)
Porcentaje de varianza acumulada
Con las primeras 11 componentes se consigue explicar aproximadamente el 90% de la varianza observada en los datos.
ggplot(data = get_eig(X = pca),
aes(x = as.factor(1:nrow(get_eig(X = pca))),
y = cumulative.variance.percent,
group = 1)
) +
geom_col(fill = "steelblue", color = "steelblue") +
geom_line(color = "black", linetype = "solid") +
geom_point(shape = 19, color = "black") +
geom_text(aes(label = paste0(round(cumulative.variance.percent, 1), "%")),
size = 3, vjust = -0.5, hjust = 0.5) +
geom_hline(yintercept = 90, color = "firebrick", linetype = "dashed") +
labs(title = "Porcentaje de varianza acumulada PCA",
x = "componente",
y = "Porcentaje de varianza acumulada") +
theme_bw()
Contribución de cada variable a las componentes
La función fviz_contrib()
del paquete factoextra
permite mostrar de forma gráfica la contribución que tiene cada variable original a una determinada componente o el promedio de contribución a un grupo de componentes.
# Contribución a la primera componente.
p1 <- fviz_contrib(X = pca, choice = "var", axes = 1, top = 10)
# Contribución a la segunda componente.
p2 <- fviz_contrib(X = pca, choice = "var", axes = 2, top = 10)
# Contribución conjunta a las 3 primeras componentes.
p3 <- fviz_contrib(X = pca, choice = "var", axes = 1:3, top = 10)
ggpubr::ggarrange(p1, p2, p3, nrow = 1)
Una vez obtenido el PCA (matriz de eigenvectors, proyecciones y medias), la reconstrucción de las observaciones iniciales se puede obtener empelando la siguiente ecuación:
\[\text{reconstrucción} = \text{PC scores} \cdot \text{Eigenvectors}^\top\] Es importante tener en cuenta que, si los datos han sido centrados y escalados (cosa que en términos generales debe hacerse al aplicar PCA) esta transformación debe revertirse.
La siguiente función obtiene la reconstrucción de las observaciones con las que se ha entrenado un objeto prcomp
, empleando únicamente determinadas componentes.
reconstruct_prcomp <- function(pca, comp = NULL){
# Esta función reconstruye las mismas observaciones con las que se ha creado el
# PCA empleando únicamente determinadas componentes principales.
# Parameters
# ----------
# pca: "prcomp"
# objeto prcomp con los resultados del PCA.
#
# comp: "numeric"
# componentes principales empleadas en la reconstrucción.
#
# Return
# ----------
# "matrix" con la reconstrucción de cada observación. Las dimensiones de la
# matriz son las mismas que las de la matriz o dataframe con el que se estrenó
# el objeto pca.
# Se comprueba que el objeto PCA es de clase "prcomp"
if (class(pca) != "prcomp") {
stop("El objeto PCA debe de haber sido creado con la función prcomp()")
}
# Si no se especifica comp, se emplean todas las componentes.
if (is.null(comp)) {
comp <- seq_along(pca$sdev)
}
# Reconstrucción
recon <- as.matrix(pca$x[, comp]) %*% t(as.matrix(pca$rotation[, comp]))
# Si se ha aplicado centrado o escalado se revierte la trasformación.
if (pca$scale[1] != FALSE) {
recon <- scale(recon , center = FALSE, scale = 1/pca$scale)
}
if (pca$center[1] != FALSE) {
recon <- scale(recon , center = -1*pca$center, scale = FALSE)
}
return(recon)
}
Se comprueba que, utilizando todas las componentes, la reconstrucción es total. Los valores reconstruidos son iguales a los datos originales (las pequeñas diferencias se deben a imprecisión numérica de las operaciones).
# Se comparan las 2 primeras variables, originales y reconstruidas, de la primera
# observación del set de datos.
print(as.numeric(reconstruccion[1, 1:2]), digits = 22)
## [1] 0.004912314661769036361338 0.693190774919389185448892
## [1] 0.004912314661768384105311 0.693190774919386298869028
El error cuadrático medio de reconstrucción de una observación se calcula como el promedio de las diferencias al cuadrado entre el valor original de sus variables y el valor reconstruido, es decir, el promedio de los errores de reconstrucción de todas sus variables elevados al cuadrado.
# Reconstrucción empleando las 11 primeras componentes (90% de la varianza explicada)
reconstruccion <- reconstruct_prcomp(pca = pca, comp = 1:11)
# Cálculo del error cuadrático medio de recostrucción
error_reconstruccion <- reconstruccion - select(datos, -y)
error_reconstruccion <- error_reconstruccion^2
error_reconstruccion <- apply(X = error_reconstruccion, MARGIN = 1, FUN = mean)
# Distribución del error de reconstrucción
tibble(error_reconstruccion) %>%
ggplot(aes(x = error_reconstruccion)) +
geom_density(fill = "steelblue") +
labs(title = "Distribución de los errores de reconstrucción (PCA)",
x = "Error de reconstrucción") +
theme_bw()
Una vez que el error de reconstrucción ha sido calculado, se puede emplear como criterio para identificar anomalías. Asumiendo que la reducción de dimensionalidad se ha realizado de forma que la mayoría de los datos (los normales) queden bien representados, aquellas observaciones con mayor error de reconstrucción deberían ser las más atípicas.
En la práctica, si se está empleando esta estrategia de detección es porque no se dispone de datos etiquetados, es decir, no se conoce qué observaciones son realmente anomalías. Sin embargo, como en este ejemplo se dispone de la clasificación real, se puede verificar si realmente los datos anómalos tienen errores de reconstrucción más elevados.
# Se añade el error de reconstrucción al dataframe original.
datos$error_reconstruccion <- error_reconstruccion
ggplot(data = datos,
aes(x = y, y = log2(error_reconstruccion))) +
geom_jitter(aes(color = y), width = 0.03, alpha = 0.3) +
geom_violin(alpha = 0) +
geom_boxplot(width = 0.2, outlier.shape = NA, alpha = 0) +
stat_summary(fun = "mean", colour = "orangered2", size = 3, geom = "point") +
labs(title = "Error de reconstruccion PCA",
x = "clasificación (0 = normal, 1 = anomalía)",
y = "log2(error reconstrucción)") +
theme_bw() +
theme(legend.position = "none")
La distribución de los errores de reconstrucción en el grupo de las anomalías (1) es claramente superior. Sin embargo, al existir solapamiento, si se clasifican las n observaciones con mayor error de reconstrucción como anomalías, se incurriría en errores de falsos positivos.
Acorde a la documentación, el set de datos Cardiotocogrpahy contiene 176 anomalías. Véase la matriz de confusión resultante si se clasifican como anomalías las 176 observaciones con mayor error de reconstrucción.
resultados <- datos %>%
select(y, error_reconstruccion) %>%
arrange(desc(error_reconstruccion)) %>%
mutate(clasificacion = if_else(row_number() <= 176, "1", "0"))
mat_confusion <- MLmetrics::ConfusionMatrix(
y_pred = resultados$clasificacion,
y_true = resultados$y
)
mat_confusion
## y_pred
## y_true 0 1
## 0 1545 110
## 1 110 66
De las 176 observaciones identificadas como anomalías, solo el 38% (66/176) lo son. El porcentaje de falsos positivos (62%) es muy elevado.
En el apartado anterior, se han empleado las 11 primeras componentes de un total de 21 ya que con ellas se puede explicar aproximadamente el 90% de la varianza observada. Que esta selección de componentes sea adecuada para la detección de anomalías depende en gran medida de en qué direcciones se desvíen las anomalías. Con frecuencia, los datos anómalos no tienen valores atípicos en las direcciones dominantes pero sí en las poco dominantes, por lo que su comportamiento atípico quedaría recogido en las últimas componentes.
Se repite la detección de anomalías pero esta vez descartando las 4 primeras componentes principales en la reconstrucción.
# Reconstrucción empleando todas excepto las 4 primeras componentes
reconstruccion <- reconstruct_prcomp(pca = pca, comp = 5:21)
# Cálculo del error cuadrático medio de reconstrucción
error_reconstruccion <- reconstruccion - select(datos, -y)
error_reconstruccion <- error_reconstruccion^2
error_reconstruccion <- apply(X = error_reconstruccion, MARGIN = 1, FUN = mean)
# Se añade el error de reconstrucción al dataframe original.
datos$error_reconstruccion <- error_reconstruccion
ggplot(data = datos,
aes(x = y, y = log2(error_reconstruccion))) +
geom_jitter(aes(color = y), width = 0.03, alpha = 0.3) +
geom_violin(alpha = 0) +
geom_boxplot(width = 0.2, outlier.shape = NA, alpha = 0) +
stat_summary(fun = "mean", colour = "orangered2", size = 3, geom = "point") +
labs(title = "Error de reconstruccion PCA",
x = "clasificación (0 = normal, 1 = anomalía)",
y = "log2(error reconstrucción)") +
theme_bw() +
theme(legend.position = "none")
resultados <- datos %>%
select(y, error_reconstruccion) %>%
arrange(desc(error_reconstruccion)) %>%
mutate(clasificacion = if_else(row_number() <= 176, "1", "0"))
mat_confusion <- MLmetrics::ConfusionMatrix(
y_pred = resultados$clasificacion,
y_true = resultados$y
)
mat_confusion
## y_pred
## y_true 0 1
## 0 1600 55
## 1 55 121
Descartando las primeras 4 componentes, 121 de las 176 (69%) observaciones identificadas como anomalías lo son realmente. Se ha conseguido reducir el porcentaje de falsos positivos al 31%.
El PCA anterior se ha obtenido empleando todas las observaciones, incluyendo potenciales anomalías. Dado que el objetivo es generar un espacio de proyección únicamente con los datos “normales”, se puede mejorar el resultado reentrenando el PCA pero excluyendo las n observaciones con mayor error de reconstrucción (potenciales anomalías).
Se repite la detección de anomalías pero esta vez descartando las observaciones con un error de reconstrucción superior al percentil 70.
# Se eliminan las observaciones con un error de reconstrucción superior al
# percentil 70%
cuantil <- quantile(x = datos$error_reconstruccion, probs = 0.7)
datos_trimmed <- datos %>% filter(error_reconstruccion < cuantil)
# Se elimina la columna error_reconstruccion de ambos dataframes para que no
# participen en el cálculo del PCA
datos$error_reconstruccion <- NULL
datos_trimmed$error_reconstruccion <- NULL
# Tras aplicar el filtrado, la columna V6 tiene varianza cero. Es necesario
# eliminarla para poder aplicar el PCA
datos$V6 <- NULL
datos_trimmed$V6 <- NULL
# PCA únicamente con los datos trimmed.
pca <- prcomp(
x = datos_trimmed %>% select(-y),
center = TRUE,
scale. = TRUE
)
Con este nuevo modelo PCA se recalcula el error de reconstrucción de todas las observaciones. La función error_reconstruccion_prcomp()
emplea un PCA ya entrenado para proyectar nuevas observaciones, reconstruirlas y calcular el error.
error_reconstruccion_prcomp <- function(pca, new_data, comp=NULL){
# Esta función calcula el error de reconstrucción de un PCA al proyectar y
# reconstruir las observaciones empleando únicamente determinadas componentes
# principales.
#
# Parameters
# ----------
# pca: "prcomp"
# objeto prcomp con los resultados del PCA
#
# new_data: "matriz" o "data.frame"
# nuevas observaciones
#
# comp: "numeric"
# componentes principales empleadas en la reconstrucción.
#
# Return
# ----------
# "numeric" vector con el error de reconstrucción de cada observación.
# Se comprueba que el objeto PCA es de clase "prcomp"
if (class(pca) != "prcomp") {
stop("El objeto PCA debe de haber sido creado con la función prcomp()")
}
# Si no se especifica comp, se emplean todas las componentes
if (is.null(comp)) {
comp <- seq_along(pca$sdev)
}
# Se seleccionan únicamente las componentes en comp
pca$rotation <- pca$rotation[, comp]
proyeciones <- predict(object = pca,
newdata = new_data
)
# Recostrucción
reconstruccion <- as.matrix(proyeciones) %*% t(as.matrix(pca$rotation))
# Si se ha aplicado centrado o escalado se revierte la trasformación
if (pca$scale[1] != FALSE) {
reconstruccion <- scale(reconstruccion , center = FALSE, scale = 1/pca$scale)
}
if (pca$center[1] != FALSE) {
reconstruccion <- scale(reconstruccion , center = -1*pca$center, scale = FALSE)
}
# Cálculo del error de reconstrucción
error_reconstruccion <- reconstruccion - new_data
error_reconstruccion <- error_reconstruccion^2
error_reconstruccion <- apply(X = error_reconstruccion, MARGIN = 1, FUN = mean)
return(error_reconstruccion)
}
# Se calcula el error de reconstrucción.
error_reconstruccion <- error_reconstruccion_prcomp(
pca = pca,
new_data = datos %>% select(-y),
comp = 5:20
)
# Se añade el error de reconstrucción al dataframe original.
datos$error_reconstruccion <- error_reconstruccion
ggplot(data = datos,
aes(x = y, y = log2(error_reconstruccion))) +
geom_jitter(aes(color = y), width = 0.03, alpha = 0.3) +
geom_violin(alpha = 0) +
geom_boxplot(width = 0.2, outlier.shape = NA, alpha = 0) +
stat_summary(fun = "mean", colour = "orangered2", size = 3, geom = "point") +
labs(title = "Error de reconstruccion PCA",
x = "clasificación (0 = normal, 1 = anomalía)",
y = "log2(error reconstrucción)") +
theme_bw() +
theme(legend.position = "none")
resultados <- datos %>%
select(y, error_reconstruccion) %>%
arrange(desc(error_reconstruccion)) %>%
mutate(clasificacion = if_else(row_number() <= 176, "1", "0"))
mat_confusion <- MLmetrics::ConfusionMatrix(
y_pred = resultados$clasificacion,
y_true = resultados$y
)
mat_confusion
## y_pred
## y_true 0 1
## 0 1605 50
## 1 50 126
Tras el reentrenamiento, 126 de las 176 (72%) observaciones identificadas como anomalías lo son realmente. Se ha reducido ligeramente el porcentaje de falsos positivos al 28%.
Al utilizar las varianzas, el método PCA es altamente sensible a outliers. Si una observación tiene un valor anómalo en una de sus variables, la varianza en esta dirección será “artificialmente” elevada. Dado que el PCA trata de identificar las direcciones con mayor varianza, el subespacio creado se habrá guiado en exceso hacia esta dirección. Existen varios métodos de PCA robusto que tratan de minimizar el impacto que tienen los datos anómalos. La mayoría de ellos sigue una de estas dos estrategias:
Obtener los eigenvectors y eigenvalues como en un PCA normal pero empleando una matriz de covarianza robusta.
Calcular directamente cada componente mediante un método robusto.
PCA basado en matriz de covarianza robusta
Esta estrategia consiste en emplear el PCA estándar pero remplazando las estimaciones clásicas de centralidad y covarianza por análogos robustos. Los dos métodos de estimación robusta más empleados son MVE (minimum volume ellipsoid) y MCD (minimum covariance determinant). Esta aproximación está implementada en la función PcaCov()
del paquete rrcov
, donde el argumento cov.control
determina el método de estimación robusta empleado.
Proyection pursuit (búsqueda de proyecciones)
Las direcciones de cada componente se buscan de forma secuencial intentando maximizar la varianza pero, en lugar de la varianza, se emplea otro estadístico de dispersión más robusto. Una ventaja de esta estrategia es que, al buscar los eigenvectors de forma secuencial, si hay muchas variables pero solo se está interesado en las primeras componentes, se ahorra tiempo de computación. Además, al no estar basado en la matriz de covarianza, se puede emplear en situaciones en las que el número de variables supera al de observaciones. Esta aproximación está implementada en la función PcaProj()
y PcaGrid()
del paquete rrcov
.
Hubert method (ROBPCA)
Esta estrategia combina las ventajas de emplear una matriz de covarianza robusta con la búsqueda de proyecciones. Su limitación principal es el tiempo de computación. Más detalles en las páginas 26-27 de An Object-Oriented Framework for Robust Multivariate Analysis, Valentin Todorov and Peter Filzmoser. Esta aproximación está implementada en la función PcaHubert()
del paquete rrcov
.
El paquete rrcov
implementa todas las estrategias de PCA robusto descritas en el apartado anterior. Se repite la detección de anomalías empleando PCA robusto con una matriz de covarianza robusta calculada con el estimador MCD.
# Cálculo de PCA
datos <- df_cardio
set.seed(123)
rpca <- PcaCov(
x = datos %>% select(-y),
cov.control = CovControlMcd(),
trace = TRUE
)
Tras obtener el PCA, se evalúa su capacidad para reducir la dimensionalidad de los datos explorando los eigenvectors, eigenvalues y la varianza capturada en el proceso. Los resultados de un PCA generado mediante PcaCov()
están almacenados de forma distinta al generado con prcomp()
, por lo que las funciones del apartado anterior no pueden emplearse.
Información sobre los eigenvalues/varianza
# Se almacenan en un tibble el valor de los eigenvalues, la varianza
# explicada y la varianza acumulada de cada componente.
rpca_eigen_var <- tibble(
componente = seq_len(rpca$k),
varianza = rrcov::getSdev(rpca)^2,
porcnt_varianza = 100 * varianza / sum(varianza),
porcnt_varianza_acumulada = cumsum(porcnt_varianza)
)
rpca_eigen_var
ggplot(data = rpca_eigen_var,
aes(x = as.factor(componente),
y = porcnt_varianza,
group = 1)
) +
geom_col(fill = "steelblue", color = "steelblue") +
geom_line(color = "black", linetype = "solid") +
geom_point(shape = 19, color = "black") +
geom_text(aes(label = paste0(round(porcnt_varianza), "%")),
size = 3, vjust = -0.5, hjust = 0.5) +
lims(y = c(0, 100)) +
labs(title = "Porcentaje de varianza por componente PCA robusto",
x = "componente",
y = "Porcentaje de varianza explicada") +
theme_bw()
Porcentaje de varianza acumulada
Con las primeras 11 componentes se consigue explicar aproximadamente el 90% de la varianza observada en los datos.
ggplot(data = rpca_eigen_var,
aes(x = as.factor(componente),
y = porcnt_varianza_acumulada,
group = 1)
) +
geom_col(fill = "steelblue", color = "steelblue") +
geom_line(color = "black", linetype = "solid") +
geom_point(shape = 19, color = "black") +
geom_text(aes(label = paste0(round(porcnt_varianza_acumulada), "%")),
size = 3, vjust = -0.5, hjust = 0.5) +
geom_hline(yintercept = 90, color = "firebrick", linetype = "dashed") +
lims(y = c(0, 100)) +
labs(title = "Porcentaje de varianza acumulada PCA robusto",
x = "componente",
y = "Porcentaje de varianza acumulada") +
theme_bw()
La siguiente función obtiene la reconstrucción de las observaciones con las que se ha entrenado un objeto PcaCov
, empleando únicamente determinadas componentes.
reconstruct_PcaCov <- function(pca, comp = NULL){
# Esta función reconstruye las mismas observaciones con las que se ha creado el
# PCA empleando únicamente determinadas componentes principales.
# Parameters
# ----------
# pca: "PcaCov"
# objeto PcaCov con los resultados del PCA.
#
# comp: "numeric"
# componentes principales empleadas en la reconstrucción.
#
# Return
# ----------
# "matrix" con la reconstrucción de cada observación. Las dimensiones de la
# matriz son las mismas que las de la matriz o dataframe con el que se estrenó
# el objeto pca.
# Se comprueba que el objeto PCA es de la clase "PcaCov"
if (class(pca) != "PcaCov") {
stop("El objeto PCA debe de haber sido creado con la función PcaCov()")
}
# Si no se especifica comp, se emplean todas las componentes.
if (is.null(comp)) {
comp <- seq_len(pca$k)
}
# Reconstrucción
recon <- as.matrix(pca$scores[, comp]) %*% t(as.matrix(pca$loadings[, comp]))
# Si se ha aplicado centrado o escalado se revierte la trasformación.
recon <- scale(recon , center = FALSE, scale = 1/pca$scale)
recon <- scale(recon , center = -1*pca$center, scale = FALSE)
return(recon)
}
Se comprueba que, utilizando todas las componentes, la reconstrucción es total. Los valores reconstruidos son iguales a los datos originales (las pequeñas diferencias se deben a imprecisión numérica de las operaciones).
# Se comparan las 2 primeras variables, originales y reconstruidas, de la primera
# observación del set de datos.
print(as.numeric(reconstruccion[1, 1:2]), digits = 22)
## [1] 0.004912314661767069184917 0.693190774919387631136658
## [1] 0.004912314661768384105311 0.693190774919386298869028
El error cuadrático medio de reconstrucción de una observación se calcula como el promedio de las diferencias al cuadrado entre el valor original de sus variables y el valor reconstruido, es decir, el promedio de los errores de reconstrucción de todas sus variables elevados al cuadrado.
# Reconstrucción empleando todas excepto las 4 primeras componentes
reconstruccion <- reconstruct_PcaCov(pca = rpca, comp = 5:20)
# Cálculo del error cuadrático medio de recostrucción
error_reconstruccion <- reconstruccion - select(datos, -y)
error_reconstruccion <- error_reconstruccion^2
error_reconstruccion <- apply(X = error_reconstruccion, MARGIN = 1, FUN = mean)
# Distribución del error de reconstrucción
tibble(error_reconstruccion) %>%
ggplot(aes(x = error_reconstruccion)) +
geom_density(fill = "steelblue") +
labs(title = "Distribución de los errores de reconstrucción PCA robusto",
x = "Error de reconstrucción") +
theme_bw()
# Se añade el error de reconstrucción al dataframe original.
datos$error_reconstruccion <- error_reconstruccion
ggplot(data = datos,
aes(x = y, y = log2(error_reconstruccion))) +
geom_jitter(aes(color = y), width = 0.03, alpha = 0.3) +
geom_violin(alpha = 0) +
geom_boxplot(width = 0.2, outlier.shape = NA, alpha = 0) +
stat_summary(fun = "mean", colour = "orangered2", size = 3, geom = "point") +
labs(title = "Error de reconstruccion PCA robusto",
x = "clasificación (0 = normal, 1 = anomalía)",
y = "log2(error reconstrucción)") +
theme_bw() +
theme(legend.position = "none")
La distribución de los errores de reconstrucción en el grupo de las anomalías (1) es claramente superior. Sin embargo, al existir solapamiento, si se clasifican las n observaciones con mayor error de reconstrucción como anomalías, se incurriría en errores de falsos positivos.
Acorde a la documentación, el set de datos Cardiotocogrpahy contiene 176 anomalías. Véase la matriz de confusión resultante si se clasifican como anomalías las 176 observaciones con mayor error de reconstrucción.
resultados <- datos %>%
select(y, error_reconstruccion) %>%
arrange(desc(error_reconstruccion)) %>%
mutate(clasificacion = if_else(row_number() <= 176, "1", "0"))
mat_confusion <- MLmetrics::ConfusionMatrix(
y_pred = resultados$clasificacion,
y_true = resultados$y
)
mat_confusion
## y_pred
## y_true 0 1
## 0 1602 53
## 1 53 123
Descartando las primeras 4 componentes y empleando PCA robusto, 123 de las 176 (70%) observaciones identificadas como anomalías lo son realmente. Se ha conseguido reducir el porcentaje de falsos positivos al 30%. Los resultados son similares a los obtenidos con PCA estándar pero sin necesidad de realizar el reentrenamiento iterativo.
Una forma de reducir la varianza de modelos estadísticos y de machine learning es el método de bagging (bootstrap aggregating). En términos generales, este método consiste entrenar n modelos, cada uno con una pseudo-muestra obtenida por resampling a partir de los datos de entrenamiento. La predicción final de una observación es el promedio de las predicciones de todos los modelos generados.
error_reconstruccion_prcomp <- function(pca, new_data, comp=NULL){
# Esta función calcula el error de reconstrucción de un PCA al proyectar y
# reconstruir las observaciones empleando únicamente determinadas componentes
# principales.
#
# Parameters
# ----------
# pca: "prcomp"
# objeto prcomp con los resultados del PCA
#
# new_data: "matriz" o "data.frame"
# nuevas observaciones
#
# comp: "numeric"
# componentes principales empleadas en la reconstrucción.
#
# Return
# ----------
# "numeric" vector con el error de reconstrucción de cada observación.
# Se comprueba que el objeto PCA es de clase "prcomp"
if (class(pca) != "prcomp") {
stop("El objeto PCA debe de haber sido creado con la función prcomp()")
}
# Si no se especifica comp, se emplean todas las componentes
if (is.null(comp)) {
comp <- seq_along(pca$sdev)
}
# Se seleccionan únicamente las componentes en comp
pca$rotation <- pca$rotation[, comp]
proyeciones <- predict(
object = pca,
newdata = new_data
)
# Recostrucción
reconstruccion <- as.matrix(proyeciones) %*% t(as.matrix(pca$rotation))
# Si se ha aplicado centrado o escalado se revierte la trasformación
if (pca$scale[1] != FALSE) {
reconstruccion <- scale(reconstruccion , center = FALSE, scale = 1/pca$scale)
}
if (pca$center[1] != FALSE) {
reconstruccion <- scale(reconstruccion , center = -1*pca$center, scale = FALSE)
}
# Cálculo del error de reconstrucción
error_reconstruccion <- reconstruccion - new_data
error_reconstruccion <- error_reconstruccion^2
error_reconstruccion <- apply(X = error_reconstruccion, MARGIN = 1, FUN = mean)
return(error_reconstruccion)
}
bagging_pca_recostruction <- function(data, iteraciones=500) {
# Esta función calcula el error de reconstrucción de las observaciones
# combinando el método de bagging con el de PCA.
# Parameters
# ----------
# data: "dataframe"
# dataframe con el valor de las variables para cada observación
#
# iteraciones: "integer"
# número de iteraciones de bagging
#
# Return
# ----------
# "vector" con el error de reconstrucción premedio de cada observación.
pb <- progress_bar$new(total = iteraciones)
mat_error <- matrix(nrow = nrow(datos), ncol = iteraciones)
for (i in 1:iteraciones) {
pb$tick()
# Resample
# --------------------------------------------------------------------------
indices_row_resample <- sample(x = 1:nrow(data), size = 0.632*nrow(data),
replace = TRUE)
indices_col_resample <- sample(x = 1:ncol(data), size = ncol(data)/3,
replace = FALSE)
resample <- data[indices_row_resample, indices_col_resample]
# Si debido al resampling, alguna columna es constante, se elimina
resample <- purrr::discard(resample, .p = function(x){var(x) == 0})
# Entrenamiento PCA
# --------------------------------------------------------------------------
pca <- prcomp(
x = resample,
center = TRUE,
scale. = TRUE
)
# Error recostrucción
# --------------------------------------------------------------------------
# Selección aleatoria de las componentes empleadas en la reconstrucción
comp_recostruccion <- sample(
x = 1:ncol(pca$rotation),
size = min(3, as.integer(ncol(pca$rotation)/3)),
replace = FALSE
)
error_reconstruccion <- error_reconstruccion_prcomp(
pca = pca,
new_data = data,
comp = comp_recostruccion
)
mat_error[, i] <- error_reconstruccion
}
# Valor promedio de todas las iteraciones
# ----------------------------------------------------------------------------
error_promedio <- apply(mat_error, 1, mean, na.rm = TRUE)
return(error_promedio)
}
datos <- df_cardio
error_reconstruccion <- bagging_pca_recostruction(
data = datos %>% select(-y),
iteraciones = 1000
)
# Distribución del error de reconstrucción
tibble(error_reconstruccion) %>%
ggplot(aes(x = error_reconstruccion)) +
geom_density(fill = "steelblue") +
labs(title = "Distribución de los errores de reconstrucción (Bagging-PCA)",
x = "Error de reconstrucción") +
theme_bw()
# Se añade el error de reconstrucción al dataframe original.
datos$error_reconstruccion <- error_reconstruccion
ggplot(data = datos,
aes(x = y, y = log2(error_reconstruccion))) +
geom_jitter(aes(color = y), width = 0.03, alpha = 0.3) +
geom_violin(alpha = 0) +
geom_boxplot(width = 0.2, outlier.shape = NA, alpha = 0) +
stat_summary(fun = "mean", colour = "orangered2", size = 3, geom = "point") +
labs(title = "Error de reconstruccion PCA",
x = "clasificación (0 = normal, 1 = anomalía)",
y = "log2(error reconstrucción)") +
theme_bw() +
theme(legend.position = "none")
resultados <- datos %>%
select(y, error_reconstruccion) %>%
arrange(desc(error_reconstruccion)) %>%
mutate(clasificacion = if_else(row_number() <= 176, "1", "0"))
mat_confusion <- MLmetrics::ConfusionMatrix(
y_pred = resultados$clasificacion,
y_true = resultados$y
)
mat_confusion
## y_pred
## y_true 0 1
## 0 1590 65
## 1 65 111
Combinando Bagging y PCA, 111 de las 176 (63%) observaciones identificadas como anomalías lo son realmente.
Los autoencoders son un tipo de redes neuronales en las que la entra y salida del modelo es la misma, es decir, redes entrenadas para predecir un resultado igual a los datos de entrada. Para conseguir este tipo de comportamiento, la arquitectura de los autoencoders es simétrica, con una región llamada encoder y otra decoder. ¿Cómo sirve esto para reducir la dimensionalidad? Los autoencoders siguen una arquitectura de cuello de botella, la región encoder está formada por una o varias capas, cada una con menos neuronas que su capa precedente, obligando así a que la información de entrada se vaya comprimiendo. En la región decoder esta compresión se revierte siguiendo la misma estructura pero esta vez de menos a más neuronas.
Para conseguir que la salida reconstruida sea lo más parecida posible a la entrada, el modelo debe aprender a capturar toda la información posible en la zona intermedia. Una vez entrenado, la salida de la capa central del autoencoder (la capa con menos neuronas) es una representación de los datos de entrada pero con una dimensionalidad igual el número de neuronas de esta capa.
La principal ventaja de los autoencoders es que no tienen ninguna restricción en cuanto al tipo de relaciones que pueden aprender, por lo tanto, a diferencia del PCA, la reducción de dimensionalidad puede incluir relaciones no lineales. La desventaja es su alto riesgo de sobreentrenamiento (overfitting), por lo que se recomienda emplear muy pocas épocas y siempre evaluar la evolución del error con un conjunto de validación.
En el caso de utilizar funciones de activación lineales, las variables generadas en el cuello de botella (la capa con menos neuronas), son muy similares a las componentes principales de un PCA pero sin que necesariamente tengan que ser ortogonales entre ellas.
El paquete h2o
permite entrenar, entre otros modelos de machine learning, redes neuronales con arquitectura de autoencoder.
# Creación de un cluster local.
h2o.init(ip = "localhost",
# Si emplea un único core los resultados son reproducibles.
# Si se dispone de muchos datos, mejor emplear varios cores.
nthreads = 1,
# Máxima memoria disponible para el cluster.
max_mem_size = "4g")
## Connection successful!
##
## R is connected to the H2O cluster:
## H2O cluster uptime: 1 hours 14 minutes
## H2O cluster timezone: Etc/UTC
## H2O data parsing timezone: UTC
## H2O cluster version: 3.32.0.1
## H2O cluster version age: 25 days
## H2O cluster name: H2O_started_from_R_ubuntu_psq559
## H2O cluster total nodes: 1
## H2O cluster total memory: 3.54 GB
## H2O cluster total cores: 8
## H2O cluster allowed cores: 1
## H2O cluster healthy: TRUE
## H2O Connection ip: localhost
## H2O Connection port: 54321
## H2O Connection proxy: NA
## H2O Internal Security: FALSE
## H2O API Extensions: Amazon S3, XGBoost, Algos, AutoML, Core V3, TargetEncoder, Core V4
## R Version: R version 4.0.3 (2020-10-10)
# Se eliminan los datos del cluster por si ya había sido iniciado.
h2o.removeAll()
h2o.no_progress()
# Se transfieren los datos al cluster de h2o.
datos <- df_cardio
datos_h2o <- as.h2o(
x = datos,
destination_frame = "datos_h2o"
)
# División de las observaciones en conjunto de entrenamiento y validación.
datos_h2o_split <- h2o.splitFrame(data = datos_h2o, ratios = 0.8, seed = 123)
datos_h2o_train <- datos_h2o_split[[1]]
datos_h2o_validacion <- datos_h2o_split[[2]]
# Se define la variables empleadas por el autoencoder
predictores <- setdiff(h2o.colnames(datos_h2o), "y")
# Entrenamiento del modelo autoencoder con 2 neuronas en la capa oculta.
autoencoder <- h2o.deeplearning(
x = predictores,
training_frame = datos_h2o_train,
validation_frame = datos_h2o_validacion,
activation = "Tanh",
autoencoder = TRUE,
hidden = c(2),
epochs = 50,
ignore_const_cols = FALSE,
score_each_iteration = TRUE,
# Seed solo válido cuando se emplea un único core
seed = 12345
)
Para identificar el número de épocas adecuado se emplea la evolución del error de entrenamiento y validación.
scoring <- as.data.frame(autoencoder@model$scoring_history)
scoring %>%
select(epochs, training_mse, validation_mse) %>%
pivot_longer(cols = c(training_mse, validation_mse),
names_to = "set",
values_to = "mse"
) %>%
ggplot() +
geom_line(aes(x = epochs, y = mse, color = set)) +
labs(title = "Evolución del error de entrenamiento y validación") +
theme_bw() +
theme(legend.position = "bottom")
A partir de las 5 épocas, la reducción en el mse es mínima. Una vez identificado el número óptimo de épocas, se reentrena el modelo, esta vez con todos los datos.
# Entrenamiento del modelo autoencoder
autoencoder <- h2o.deeplearning(
x = predictores,
training_frame = datos_h2o,
activation = "Tanh",
autoencoder = TRUE,
hidden = c(2),
epochs = 5,
ignore_const_cols = FALSE,
score_each_iteration = TRUE,
# Seed solo válido cuando se emplea un único core
seed = 12345
)
La función h2o.anomaly()
permite obtener el error de reconstrucción empleando un modelo H2OAutoEncoderModel
. Para ello, realiza automáticamente la codificación, decodificación y la comparación de los valores reconstruidos con los valores originales.
error_reconstruccion <- h2o.anomaly(
object = autoencoder,
data = datos_h2o[, predictores],
per_feature = FALSE
)
error_reconstruccion <- as.data.frame(error_reconstruccion)
error_reconstruccion <- as.data.frame(error_reconstruccion)
error_reconstruccion <- error_reconstruccion[, 1]
# Distribución del error de reconstrucción.
tibble(error_reconstruccion) %>%
ggplot(aes(x = error_reconstruccion)) +
geom_density(fill = "steelblue") +
labs(title = "Distribución de los errores de reconstrucción (Autoencoder)",
x = "Error de reconstrucción") +
theme_bw()
Una vez que el error de reconstrucción ha sido calculado, se puede emplear como criterio para identificar anomalías. Asumiendo que la reducción de dimensionalidad se ha realizado de forma que la mayoría de los datos (los normales) queden bien representados, aquellas observaciones con mayor error de reconstrucción deberían ser las más atípicas.
En la práctica, si se está empleando esta estrategia de detección es porque no se dispone de datos etiquetados, es decir, no se conoce qué observaciones son realmente anomalías. Sin embargo, como en este ejemplo se dispone de la clasificación real, se puede verificar si realmente los datos anómalos tienen errores de reconstrucción más elevados.
# Se añade el error de reconstrucción al dataframe original.
datos$error_reconstruccion <- error_reconstruccion
ggplot(data = datos,
aes(x = y, y = log2(error_reconstruccion))) +
geom_jitter(aes(color = y), width = 0.03, alpha = 0.3) +
geom_violin(alpha = 0) +
geom_boxplot(width = 0.2, outlier.shape = NA, alpha = 0) +
stat_summary(fun = "mean", colour = "orangered2", size = 3, geom = "point") +
labs(title = "Error de reconstruccion Autoencoder",
x = "clasificación (0 = normal, 1 = anomalía)",
y = "log2(error reconstrucción)") +
theme_bw() +
theme(legend.position = "none")
Acorde a la documentación, el set de datos Cardiotocogrpahy contiene 176 anomalías. Véase la matriz de confusión resultante si se clasifican como anomalías las 176 observaciones con mayor error de reconstrucción.
resultados <- datos %>%
select(y, error_reconstruccion) %>%
arrange(desc(error_reconstruccion)) %>%
mutate(clasificacion = if_else(row_number() <= 176, "1", "0"))
mat_confusion <- MLmetrics::ConfusionMatrix(
y_pred = resultados$clasificacion,
y_true = resultados$y
)
mat_confusion
## y_pred
## y_true 0 1
## 0 1571 84
## 1 84 92
De las 176 observaciones identificadas como anomalías, solo el 52% (92/176) lo son realmente. El porcentaje de falsos positivos es del 48%.
El autoencoder anterior se ha entrenado empleando todas las observaciones, incluyendo las potenciales anomalías. Dado que el objetivo es generar un espacio de proyección para datos “normales”, se puede mejorar el resultado reentrenando el modelo pero esta vez excluyendo las n observaciones con mayor error de reconstrucción (potenciales anomalías).
# Se eliminan las observaciones con un error de reconstrucción superior al
# percentil 70
cuantil <- quantile(x = error_reconstruccion, probs = 0.7)
datos_trimmed <- datos %>% filter(error_reconstruccion < cuantil)
autoencoder <- h2o.deeplearning(
x = predictores,
training_frame = datos_trimmed_h2o,
activation = "Tanh",
autoencoder = TRUE,
hidden = 2,
epochs = 5,
score_each_iteration = TRUE,
ignore_const_cols = FALSE,
# Seed solo válido cuando se emplea un único core
seed = 12345
)
error_reconstruccion <- h2o.anomaly(
object = autoencoder,
data = datos_h2o[, predictores],
per_feature = FALSE
)
error_reconstruccion <- as.data.frame(error_reconstruccion)
error_reconstruccion <- as.data.frame(error_reconstruccion)
error_reconstruccion <- error_reconstruccion[, 1]
# Se añade el error de reconstrucción al dataframe original.
datos$error_reconstruccion <- error_reconstruccion
ggplot(data = datos,
aes(x = y, y = log2(error_reconstruccion))) +
geom_jitter(aes(color = y), width = 0.03, alpha = 0.3) +
geom_violin(alpha = 0) +
geom_boxplot(width = 0.2, outlier.shape = NA, alpha = 0) +
stat_summary(fun = "mean", colour = "orangered2", size = 3, geom = "point") +
labs(title = "Error de reconstruccion Autoencoder",
x = "clasificación (0 = normal, 1 = anomalía)",
y = "log2(error reconstrucción)") +
theme_bw() +
theme(legend.position = "none")
resultados <- datos %>%
select(y, error_reconstruccion) %>%
arrange(desc(error_reconstruccion)) %>%
mutate(clasificacion = if_else(row_number() <= 176, "1", "0"))
mat_confusion <- MLmetrics::ConfusionMatrix(
y_pred = resultados$clasificacion,
y_true = resultados$y
)
mat_confusion
## y_pred
## y_true 0 1
## 0 1590 65
## 1 65 111
Tras el reentrenamiento, 111 de las 176 (63%) observaciones identificadas como anomalías lo son realmente. Se ha conseguido reducir el porcentaje de falsos positivos al 37%.
sesion_info <- devtools::session_info()
dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
)
Outlier Analysis Aggarwal, Charu C.
Deep Learning with H2O by Arno Candel & Erin LeDell with assistance from Viraj Parmar & Anisha Arora Edited by: Angela Bartz
¿Cómo citar este documento?
Detección de anomalías: Autoencoders y PCA por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/52_deteccion_anomalias_autoencoder_pca.html
¿Te ha gustado el artículo? Tu ayuda es importante
Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊
Este material, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.
Se permite:
Compartir: copiar y redistribuir el material en cualquier medio o formato.
Adaptar: remezclar, transformar y crear a partir del material.
Bajo los siguientes términos:
Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.
NoComercial: No puedes utilizar el material para fines comerciales.
CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.