Resultados

Analisis exploratorio

Pequeño y feo informe en RPubs.

Preprocesado

Fueron suministrados dos conjuntos de datos:

Datos dia 1: 40 espectros por muestra. (40 * 5 = 200 espectros)
Datos dia 2: 100 espectros por muestra. ( 100 * 5 = 500 espectros)

Los conjuntos se mantendrán separados (la razón se explica mas adelante) y se realizan los siguientes pasos de pre procesado:

Promediar 3 espectros: Con esto se intenta disminuir la dispersión en las observaciones, sin embargo, hay que tener en cuenta que el numero de observaciones se reduce a 1/3 del tamaño original, pasando de 40 a 13, y de 100 a 33, para cada muestra.
Según lo observado en el análisis exploratorio, solo se usaron datos del detector numero 2.
Normalizar espectros por suma total: Referencia: Impact of Laser-Induced Breakdown Spectroscopy data normalization on multivariate classification accuracy

# librerias usadas
library(tidyverse)      # para manipulaciond de los datos
#setwd('./deteccion de H')
source('./funciones utiles.R')
source('./BaseLine_Script.R')
source('./Funciones_Prepocesado.R')

data_path_1 <- "./Data/Calibracion Zr2.5Nb - 4.53 J- 2.92us/"
data_path_2 <- "./Data/new data/"

carpetas <- list('ARG-2','ARG-4','ARG-3','ARG-5','ARG-6')

# Matriz con 40 espectros por muestra
listaM.1 <- lapply(carpetas, function(x){df_func(data_path_1, x, 40)})
listaM.2 <- lapply(carpetas, function(x){df_func(data_path_2, x, 100)})

# datos correspondientes al detector 2
listaM.1 <- map(listaM.1, ~ .x %>% .[,2049:3983] %>%  data.frame() )
listaM.2 <- map(listaM.2, ~ .x %>% .[,2049:3983] %>%  data.frame() )

names(listaM.1) <- c("2", "23", "42","79","99")
names(listaM.2) <- c("2", "23", "42","79","99")

rm(carpetas, data_path_1, data_path_2)

# Función para reordenar los datos de forma aleatoria
FUN.random <- function(M){
        size <- nrow(M)
        index <- sample(seq_len(size), size = size)
        M <- M[index,]
        M }

listaM.1 <- map(listaM.1, FUN.random)
listaM.2 <- map(listaM.2, FUN.random)

# Promediar 3 espectros, ahora tenemos 33 y 13
listaM.1 <- map(listaM.1, FUN.prom.4spec, i=3)
listaM.2 <- map(listaM.2, FUN.prom.4spec, i=3)

# Normalizacion por suma total 
listaM.1 <- map(listaM.1, ~ apply(.x, 1, FUN.norm.spec) %>% t() )
listaM.2 <- map(listaM.2, ~ apply(.x, 1, FUN.norm.spec) %>% t() )

Finalmente se sumaron 3 longitudes de onda con el objetivo de disminuir la cantidad de variables que describen a cada observación. Se disminuyen las características a 1/3 del numero original, esto es importante dado que cada observación (espectro) tiene un numero elevado de características si se compara con el numero de observaciones disponibles.

listaM.1 <- map(listaM.1, ~ apply(.x, 1, FUN.sum.long, n = ncol(.x)/3 ) %>% t() )
listaM.2 <- map(listaM.2, ~ apply(.x, 1, FUN.sum.long, n = ncol(.x)/3 ) %>% t() )

Sustracción de la linea base: Referencias: Influence of baseline subtraction on laser-induced breakdown spectroscopic data. Metodo utilizado para hallar linea base

L.1 <- map(listaM.1, ~ .x %>% apply(1, BaseLine, w= 25) )
L.1[[1]] %>% plot.comparison() #%>% plotly::ggplotly()

L.2 <- map(listaM.2, ~ .x %>% apply(1, BaseLine, w= 25) )
L.2[[1]] %>% plot.comparison() #%>% plotly::ggplotly()

El objetivo general de este modelado, como en cualquier otro problema de Machine Learning, es la generalización; en el mejor de los casos un modelo entrenado con los datos del día 2 (dado que es conjunto de mayor tamaño), aproximaría correctamente los datos del día 1, tanto en un problema de regresión como de clasificación. Por esta razón se mantendrán separados los datos del día 1 y día 2.

Todo el código del preprocesado (y los modelos) esta en mi repositorio de GitHub

Creación de predictores

Se ajusto visualmente un limite para la selección de longitudes de onda de interés.

I_raw <- L.1 %>% 
        # seleccionar de c/muestra y c/espectro  la columna 'I' 
        map(~.x %>% map_dfc(~.x %>% select('I')) %>% t(), .x ) %>% 
        # Promediar para cada muestra
        map(~ apply(.x, 2, mean)) %>% 
        # Todo a un solo dataframe
        bind_rows()

BL <- L.1 %>% 
        # seleccionar de c/muestra y c/espectro  la columna 'Int.corrected' 
        map(~.x %>% map_dfc(~.x %>% select('Bi')) %>% t(), .x ) %>% 
        # Promediar para cada muestra
        map(~ apply(.x, 2, mean)) %>% 
        # Todo a un solo dataframe
        bind_rows()

# Ajuste visual de Z
z <- 1.25

g <- tibble(I = I_raw$`99`, Base = BL$`99`*z) %>% 
        rowid_to_column() %>% 
        ggplot(aes(rowid, I)) + 
        geom_line() +
        geom_line(aes(rowid, Base), colour = 'red') 
g

Gráficamente, todos los puntos que sobrepasen la linea roja son las nuevas características que describen el espectro. (Seria lo mismo pasar una linea recta sobre el espectro ya corregido, pero en fases previas usé este código para ver si era necesario sustraer la linea base, así que quedo así, luego se corregirá)

*OJO el desplazamiento de esta linea es un posible parámetro a tunear en futuros modelos

Este procedimiento se aplica para un espectro promedio de cada muestra, 2ppm, 23ppm, 42ppm, 79ppm, 99ppm; con esto se generan características que lo describen y que no son necesariamente las mismas para cada muestra:

BL.z <- BL * z
indices <- I_raw > BL.z
colnames(indices) <- c('a','b','c','d','f')
indices <- indices %>% as.data.frame() %>% mutate(index = ifelse((a == TRUE) | (b  == TRUE) | (c  == TRUE) | (d  == TRUE) | (f  == TRUE), TRUE, FALSE))

Para la muestra de 2ppm se obtienen asi 212, para la muestra de 23ppm se obtienen 210, para la muestra de 42ppm se obtienen 211, para la muestra de 79ppm se obtienen 214 y para la muestra de 99ppm se obtienen 247. Finalmente, al unir los indices que describen a cada muestra se obtienen un total de 247 caracteristicas de interes (predictores).

Finalmente, seleccionamos las caracteristicas de interés:

# Seleccionando características de interés
I_new.1 <- L.1 %>% 
        # seleccionar de c/muestra y c/espectro  la columna 'Int.corrected' 
        map(~.x %>% map_dfc(~.x %>% select('Int.corrected')) %>% t(), .x ) %>% 
        map(~.x %>% as.data.frame()) %>% 
        bind_rows(.id = 'class')

df_temp <- I_new.1 %>% select(V1:V645)
df_temp <- df_temp[,indices$index]

data.1 <- data.frame(class = I_new.1$class, df_temp ) # 13 * 5 = 65 rows

I_new.2 <- L.2 %>% 
        # seleccionar de c/muestra y c/espectro  la columna 'Int.corrected' 
        map(~.x %>% map_dfc(~.x %>% select('Int.corrected')) %>% t(), .x ) %>% 
        map(~.x %>% as.data.frame()) %>% 
        bind_rows(.id = 'class')

df_temp <- I_new.2 %>% select(V1:V645)
df_temp <- df_temp[,indices$index]

data.2 <- data.frame(class = I_new.2$class, df_temp )

Para el dia 1 se cuenta con 65 observaciones y su distribucion es:

## 
##  2 23 42 79 99 
## 13 13 13 13 13

Para el dia 2 se cuenta con 165 observaciones y su distribucion es:

## 
##  2 23 42 79 99 
## 33 33 33 33 33

Estos son los datos que posteriormente se distribuirán en los tres conjuntos llamados train, val y test.

Modelado

El modelado se centra en el uso de redes neuronales (de tipo Feed Forward Neural Network), pero se contempla experimentar con otros métodos, así como distintos tipos de pre procesados (en curso). Se presentan 3 modelos que comparten los mismos pasos de pre procesamiento de los datos antes mencionados.

Es practica común dividir los datos en en dos grupos un “train set” y un “test set” (80% y 20% respectivamente). Para aumentar la probabilidad de generalizar correctamente, dada.. bueno… la confiabilidad en la situación del LIBS, los datos se dividieron en 3 sets; train (74%), validation (13%) y test (13%). El modelo se entreno usando el test set y se optimizo internamiente utilizando Cross Validation. Sin embargo, los modelos obtenidos se eligieron en función a su desempeño con el validation set y no a su desempeño en el paso de Cross Validation ya que se corre el riesgo de memorizar los datos de entrada y alejarnos de la generalización.

Para la optimizacion de los hiperparametros se utilizo random grid search.

Modelo 1

Este modelo plantea el escenario ideal en el que se toman datos en un día cualquiera, se entrena un modelo, y luego se evalúa en un nuevo set de datos. Se plantea un modelo de clasificación porque el entrenamiento es mas rápido en mi limitada compu :(. La intensión es observar como afecta la distribución de los datos al comportamiento final, luego se pasara a un modelo de regresión.

La distribución de los datos fue la siguiente:

Datos dia 1: Total (por muestra) = 13 -> Train = 0, validation = 0, test = 13
Datos dia 2: Total (por muestra) = 33 -> Train = 25, validation = 4, test = 4

El codigo esta en mi repo, acá se carga directamente el mejor modelo y se evalúa su desempeño en el test set.

load(file = './Data/Data_modelo_C.RData')

library(h2o)    # Libreria de machine learning        
h2o.init(nthreads = 2)

##  Connection successful!
## 
## R is connected to the H2O cluster: 
##     H2O cluster uptime:         3 hours 27 minutes 
##     H2O cluster timezone:       America/Argentina/Buenos_Aires 
##     H2O data parsing timezone:  UTC 
##     H2O cluster version:        3.32.1.7 
##     H2O cluster version age:    1 month and 4 days  
##     H2O cluster name:           H2O_started_from_R_gomez_hcv607 
##     H2O cluster total nodes:    1 
##     H2O cluster total memory:   3.87 GB 
##     H2O cluster total cores:    4 
##     H2O cluster allowed cores:  2 
##     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, Algos, AutoML, Core V3, TargetEncoder, Core V4 
##     R Version:                  R version 4.1.1 (2021-08-10)

h2o.no_progress()

train <- as.h2o(train)
val <- as.h2o(val)
test <- as.h2o(test)

train[,'class'] <- as.factor(train[,'class'])
val[,'class'] <- as.factor(val[,'class'])
test[,'class'] <- as.factor(test[,'class'])

mdl1 <- h2o.loadModel("./outputs/model_C/DeepLearning_grid__2_AutoML_20211006_110049_model_29")
h2o.confusionMatrix(mdl1)

## Confusion Matrix: Row labels: Actual class; Column labels: Predicted class
##         2 23 42 79 99  Error      Rate
## 2      25  0  0  0  0 0.0000 =  0 / 25
## 23      0 25  0  0  0 0.0000 =  0 / 25
## 42      0  0 25  0  0 0.0000 =  0 / 25
## 79      0  0  0 25  0 0.0000 =  0 / 25
## 99      0  0  0  0 25 0.0000 =  0 / 25
## Totals 25 25 25 25 25 0.0000 = 0 / 125

h2o.confusionMatrix(mdl1, val)

## Confusion Matrix: Row labels: Actual class; Column labels: Predicted class
##        2 23 42 79 99  Error     Rate
## 2      4  0  0  0  0 0.0000 =  0 / 4
## 23     0  4  0  0  0 0.0000 =  0 / 4
## 42     0  0  4  0  0 0.0000 =  0 / 4
## 79     0  0  0  4  0 0.0000 =  0 / 4
## 99     0  0  0  0  4 0.0000 =  0 / 4
## Totals 4  4  4  4  4 0.0000 = 0 / 20

h2o.confusionMatrix(mdl1, test)

## Confusion Matrix: Row labels: Actual class; Column labels: Predicted class
##        2 23 42 79 99  Error      Rate
## 2      6  6  0  0  5 0.6471 = 11 / 17
## 23     0  4  0  3 10 0.7647 = 13 / 17
## 42     0  8  5  2  2 0.7059 = 12 / 17
## 79     0  0  8  7  2 0.5882 = 10 / 17
## 99     0  0  1  0 16 0.0588 =  1 / 17
## Totals 6 18 14 12 35 0.5529 = 47 / 85

Claramente este modelo no puede generalizar nada bien, pero, ahora separo los datos de acuerdo al conjunto al cual pertenecen (dia 1 o dia 2):

real.ppm <- as.data.frame(test[,'class']) 
mdl1.ppm <- as.data.frame(h2o.predict(mdl1, test))
df <- data.frame(real.ppm, mdl1.ppm[,'predict']) %>% 
        setNames(c('real','mdl1')) %>% 
        mutate( error1 = real == mdl1 )

# Que pasa si separo el test set por dia? cual es el error para cada dia?
Resul.d1 <- df[1:65,]
Resul.d2 <- df[66:85,]

mean.dia1 <- mean(Resul.d1$real  ==  Resul.d1$mdl1)     # 0.6
mean.dia2 <- mean(Resul.d2$real  ==  Resul.d2$mdl1)     # 0.95

El porcentaje de aciertos para el dia 1 es mean.dia1, y se distribuyen de la siguiente manera:

# separando por dia y por muestra
Resul.d1[,c('real','mdl1')] %>% 
        mutate(error = real == mdl1) %>%
        group_by(real) %>% 
        summarise(mean = mean(error))

## # A tibble: 5 x 2
##   real    mean
##   <fct>  <dbl>
## 1 2     0.154 
## 2 23    0     
## 3 42    0.0769
## 4 79    0.308 
## 5 99    0.923

El porcentaje de aciertos para el dia 1 es mean.dia2, y se distribuyen de la siguiente manera:

Resul.d2[,c('real','mdl1')] %>% 
        mutate(error = real == mdl1) %>%
        group_by(real) %>% 
        summarise(mean = mean(error))

## # A tibble: 5 x 2
##   real   mean
##   <fct> <dbl>
## 1 2      1   
## 2 23     1   
## 3 42     1   
## 4 79     0.75
## 5 99     1

El modelo no es bueno generalizando. ¡Uh! que dolor, que dolor, que pena.

rm(list= ls())

Modelo 2

Se cambia la distribución de los datos:

Datos dia 1: Total (por muestra) = 13 -> Train = 9, validation = 2, test = 2
Datos dia 2: Total (por muestra) = 33 -> Train = 25, validation = 4, test = 4

El código esta en mi repo, acá se cargan directamente los mejores modelos y se evalúa su desempeño en el test set.

Modelo 2 - clasificacion:

load(file = './Data/Data_modelo_A.RData')

# library(h2o)    # Libreria de machine learning        
# h2o.init(nthreads = 2)

train <- as.h2o(train)
val <- as.h2o(val)
test <- as.h2o(test)

train[,'class'] <- as.factor(train[,'class'])
val[,'class'] <- as.factor(val[,'class'])
test[,'class'] <- as.factor(test[,'class'])

mdl1 <- h2o.loadModel("./outputs/DeepLearning_grid__1_AutoML_20211004_115311_model_31")
h2o.confusionMatrix(mdl1)

## Confusion Matrix: Row labels: Actual class; Column labels: Predicted class
##         2 23 42 79 99  Error      Rate
## 2      34  0  0  0  0 0.0000 =  0 / 34
## 23      0 34  0  0  0 0.0000 =  0 / 34
## 42      0  0 34  0  0 0.0000 =  0 / 34
## 79      0  0  0 34  0 0.0000 =  0 / 34
## 99      0  0  0  0 34 0.0000 =  0 / 34
## Totals 34 34 34 34 34 0.0000 = 0 / 170

h2o.confusionMatrix(mdl1, val)

## Confusion Matrix: Row labels: Actual class; Column labels: Predicted class
##        2 23 42 79 99  Error     Rate
## 2      6  0  0  0  0 0.0000 =  0 / 6
## 23     0  6  0  0  0 0.0000 =  0 / 6
## 42     0  0  6  0  0 0.0000 =  0 / 6
## 79     0  0  0  6  0 0.0000 =  0 / 6
## 99     0  0  0  0  6 0.0000 =  0 / 6
## Totals 6  6  6  6  6 0.0000 = 0 / 30

h2o.confusionMatrix(mdl1, test)

## Confusion Matrix: Row labels: Actual class; Column labels: Predicted class
##        2 23 42 79 99  Error     Rate
## 2      6  0  0  0  0 0.0000 =  0 / 6
## 23     0  4  0  2  0 0.3333 =  2 / 6
## 42     0  0  6  0  0 0.0000 =  0 / 6
## 79     0  0  0  5  1 0.1667 =  1 / 6
## 99     0  0  0  0  6 0.0000 =  0 / 6
## Totals 6  4  6  7  7 0.1000 = 3 / 30

El modelo ahora tiene 90% de aciertos en el test set, pero, eso no significa que pueda generalizar.

OJO los datos del validation set soy muy representativos del training set, y como se vio en la matriz de confusión anterior, también lo son del test set. Pero ¿que pasaría si tomamos datos un tercer día, la variación de los datos interdia (que no podemos medir, los famosos desconocidos desconocidos, o conocidos desconocidos como energía del láser, humedad, distancia) se ve lo suficientemente representada por dos días de muestreo?

rm(list= ls())

Modelo 2 - regresion:

Misma distribución de datos, pero ahora como problema de regresión.

load(file = './Data/Data_modelo_A.RData')
# library(h2o)    # Libreria de machine learning        
# h2o.init(nthreads = 2)

train <- as.h2o(train)
val <- as.h2o(val)
test <- as.h2o(test)

train$class <- as.numeric(train$class) 
val$class <- as.numeric(val$class) 
test$class <- as.numeric(test$class) 

mdl1 <- h2o.loadModel("C:\\Users\\gomez\\Documents\\LIBS\\deteccion de H\\outputs\\DeepLearning_grid__3_AutoML_20211004_154141_model_8")
mdl2 <- h2o.loadModel("C:\\Users\\gomez\\Documents\\LIBS\\deteccion de H\\outputs\\DeepLearning_grid__3_AutoML_20211004_154141_model_3")
mdl3 <- h2o.loadModel("C:\\Users\\gomez\\Documents\\LIBS\\deteccion de H\\outputs\\DeepLearning_grid__2_AutoML_20211004_154141_model_8")

Se guardaron los tres mejores modelos para futura comparación con nuevos datos.

real.ppm <- as.data.frame(test[,'class']) 
mdl1.ppm <- as.data.frame(h2o.predict(mdl1, test))
mdl2.ppm <- as.data.frame(h2o.predict(mdl2, test))
mdl3.ppm <- as.data.frame(h2o.predict(mdl3, test))

error <- data.frame(real.ppm, mdl1.ppm, mdl2.ppm, mdl3.ppm) %>% 
        summarise(diff1 = abs(class - predict),
                  diff2 = abs(class - predict.1),
                  diff3 = abs(class - predict.2))

Analisando los errores absolutos para cada modelo:

error %>% map_df(quantile, c(0.25,0.5,0.75,0.95, 1))

## # A tibble: 3 x 5
##   `25%` `50%` `75%` `95%` `100%`
##   <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 0.687 0.805  1.31 14.2    19.4
## 2 0.268 0.951  4.01 14.1    30.1
## 3 0.439 0.930  2.86  9.46   32.0

OJO Son modelos interesantes para probar con datos de un tercer día.

Modelo 3

Finalmente, se cambia la distribución de los datos para un tercer modelado:

Datos dia 1: Total (por muestra) = 13 -> Train = 0, validation = 6, test = 7
Datos dia 2: Total (por muestra) = 33 -> Train = 25, validation = 0, test = 8

El objetivo nuevamente es la generalización. Esta vez el modelo se entrena solo con datos del día 2 pero se evalúa exclusivamente con datos del día 1, por lo tanto no aprende directamente de los datos del día 1. Finalmente se evalúan tres modelos contra un test set mixto dia 1 + dia 2.

El código esta en mi repo, acá se cargan directamente los mejores modelos y se evalúa su desempeño en el test set.

Modelo 3 - clasificacion:

train <- as.h2o(train)
val <- as.h2o(val)
test <- as.h2o(test)

train[,'class'] <- as.factor(train[,'class'])
val[,'class'] <- as.factor(val[,'class'])
test[,'class'] <- as.factor(test[,'class'])

mdl1 <- h2o.loadModel("C:\\Users\\gomez\\Documents\\LIBS\\deteccion de H\\outputs\\DeepLearning_grid__3_AutoML_20211005_201843_model_3")

predictions <- h2o.predict(mdl1, test)
real.class <- as.data.frame(test[,'class']) 
pred.class <- as.data.frame(predictions)
df <- data.frame(real.class,pred.class[,'predict'])
names(df) <- c('real','predicted')

mean.all <- mean((as.numeric(real.class$class) == as.numeric(pred.class$predict))) # 0.78!!!

Este modelo tiene un porcentaje de aciertos de mean.all en el test set!! Pero, que pasa si evaluamos por dia?

# Que pasa si separo el test set por dia? cual es el error para cada dia?
Resul.d1 <- df[1:35,]
Resul.d2 <- df[36:75,]

mean.d1 <- mean(Resul.d1$real  ==  Resul.d1$predicted)     # 0.6
mean.d2 <- mean(Resul.d2$real  ==  Resul.d2$predicted)     # 0.95

El porcentaje de aciertos para los datos del test set provenientes del dia 1 es mean.d1, separados por muestra obtenemos:

Resul.d1 %>% 
        mutate(error = real == predicted) %>%
        group_by(real) %>% 
        summarise(mean = mean(error))

## # A tibble: 5 x 2
##   real   mean
##   <fct> <dbl>
## 1 2     0.571
## 2 23    0.286
## 3 42    0.429
## 4 79    1    
## 5 99    0.714

El porcentaje de aciertos para los datos del test set provenientes del dia 2 es mean.d2, separados por muestra obtenemos:

Resul.d2 %>% 
        mutate(error = real == predicted) %>%
        group_by(real) %>% 
        summarise(mean = mean(error))

## # A tibble: 5 x 2
##   real   mean
##   <fct> <dbl>
## 1 2      1   
## 2 23     1   
## 3 42     1   
## 4 79     0.75
## 5 99     1

Acá se ve algo muy interesante…

Modelo 3 - regresion: