Tyto alba (brnowl)

Búho pálido, blanco abajo y anaranjado con manchas grises arriba. Disco facial en forma de corazón blanco. Caza principalmente roedores en áreas abiertas por la noche. Duerme en edificios viejos y cajas de nidos durante el día. Ampliamente distribuido en todos los continentes, excepto en la Antártica. Por la noche, parece muy blanco en faros. Su llamado es apropiado para Halloween: un chillido escalofriante.

Actitis hypoleucos (comsan)

Migradora parcial en la península ibérica con individuos invernantes, en paso e incluso veraneantes procedentes de las poblaciones europeas del norte y centro. Escasa población reproductora. Aparecen en una amplia variedad de hábitat acuáticos, pero selecciona para reproducirse orillas pedregosas de ríos o lagos. Generalmente individuos solitarios, pero también pueden observarse pequeños grupos. Evita bandos mixtos con otras especies. Camina con un distintivo y constante movimiento “arriba-abajo” de la cola, tipo wagtail. Tiene plumaje pardusco en el dorso, donde destaca el plumaje blanco del hombro que se ve como una pequeña franja blanca entre el ala plegada y el pecho.

Anas platyrhynchos (mallar3)

El típico pato dentro de su rango, se encuentra en cualquier lugar donde hay agua, incluyendo parques de ciudades, arroyos en el jardín y diferentes humedales. Los machos tienen cabeza verde, pecho castaño y cuerpo gris. Las hembras son marrones moteadas con manchas anaranjadas y negras en el pico. Las barras blancas enfrente y atrás del parche azul del ala son más conspicuas comparadas a las de American Black Duck y el Mottled Duck.

Alauda arvensis (skylar)

Bastante común en campos abiertos, especialmente pastizales en terrenos escabrosos, páramos, tierras de cultivo, aeropuertos. Como la mayoría de los aláudidos, normalmente es poco visible en el suelo y mejor detectada por la voz. Posee un canto prolongado en vuelo, compuesto de gorjeos y trinos, a menudo tan alto que apenas se ve una mancha en el cielo, si es que está visible. Si está posada, destaca el plumaje rayado marrón, cresta tupida distintiva, y anillo ocular ancho, difuso y pálido. Más corpulenta que un bisbita pero y a diferencia de ésta no suele menear su parte trasera. En vuelo muestra lados blancos de la cola, así como un borde posterior blanco y estrecho en las alas.

Cardinalis cardinalis (norcar)

El cardenal norteño es un pájaro cantor de tamaño medio con una longitud corporal de 21-23 cm. Tiene un distintivo penacho y una máscara en la cara que es negra en el macho y gris en la hembra. Presenta dimorfismo sexual en su coloración; el plumaje del macho es rojo brillante, mientras que el de la hembra es de un tono opaco que mezcla rojo y café. El cardenal norteño es predominantemente granívoro, pero también se alimenta de insectos y fruta. El macho tiene un comportamiento territorial, delimitando su territorio con su canto

Passer domesticus (houspa)

Ampliamente distribuido y abundante en ciudades, vecindarios y granjas. Evita bosques densos. Las bandadas se agrupan en densos arbustos, moviéndose y vocalizando entre ellos. Los machos tienen baberos negros elegantes, nuca rufa brillante, y alas con brillantes tonos marrones y caqui. Las partes inferiores son gris pálido. Las hembras son marrón uniforme con una cara linda y ceja más clara. Nativo de Eurasia; introducido a Norteamérica y no está estrechamente relacionado con otros gorriones de allí.

Filtramos la data para obtener los registros de las especies que nos interesan, como también aquellas columnas con información relevante:

sub_data <- data %>% 
  select(-c('author', 'license','url',
                                 'filename')) %>% 
  filter(primary_label %in% birds)
head(sub_data) %>% kable()

# c('secondary_labels','type', 'author', 'license','url',
#                                  'filename')

sub_data <- sub_data %>% 
  rename(etiqueta=primary_label, 
         etiqueta2=secondary_labels,
         latitud=latitude,
         type=type,
         longitud=longitude,
         Nombre_cientifico=scientific_name,
         Nombre_comun=common_name,
         tiempo=time)
head(sub_data) %>%  kable()

Se tienen 3000 registros en total de seis especies de aves:

dim(sub_data)

## [1] 3000    9

Primary_label

Inicialmente se tienen registros de espectrogramas para 152 especies de aves. La distribución de datos para cada especie se puede observar a continuación:

# Frecuencia de cada especie de ave
library(forcats)
fig <- data %>% 
  count(primary_label) %>% 
  mutate(primary_label = fct_reorder(primary_label, n, .desc = T)) %>% 
    plot_ly(x=~primary_label, y=~n,
            marker = list(color = 'rgb(118, 215, 196)',
                      line = list(color = 'rgb(14, 102, 85)',
                                  width = 1.5))) %>% 
  add_bars()
fig %>% 
  layout(title = list(text="Frecuencia de Especies de aves", y = 0.95),
         xaxis = list(title = 'Especies'),
         yaxis = list(title = 'Frecuencia'))

Solo hay 6 especies de aves con 500 registros. Por esta razón, y la imposibilidad actual de aplicar técnicas de aumetación de datos para los datos, es que se decide trabajar con estás 6 especies , ya que son las que tienen mayor cantidad de datos.

De aquí en adelante los análisis se realizarán solo para estás 6 especies.

secondary_label

La etiqueta secundaria es una etiqueta asignada por la persona que generó la grabación, indicando que el audio podría tener más de una etiqueta o clasificación, debido a la presencia de otras sonidos de aves en la grabación, principalmente aves que se encuentran en el trasfondo de la grabación.

# Determinar si existe una etiqueta secundaria
sub_data$etiqueta2 <- str_replace_all(sub_data$etiqueta2, "[\\[\\]']", "")
sub_data$et2_count <- lengths(strsplit(sub_data$etiqueta2, ','))
sub_data %>% 
  group_by(et2_count) %>% 
  summarise(conteo = n()) %>%  kable()

La mayoría de audios no registran sonidos de más de una especie de ave.

rating

La variable de rating permite identificar la calidad que tiene el audio con referencia al sonido emitido por el ave. En el siguiente grafico se puede observar, que la mayoría de audios tienen un puntaje de 4, seguido por un puntaje de 5, lo que indica que los audios con los que se trabajará son adecuados para entrenar el modelo predictivo.

library(forcats)

fig <- sub_data %>% 
  count(rating) %>% 
  # mutate(rating = fct_reorder(rating, n, .desc = T)) %>% 
    plot_ly(x=~rating, y=~n, 
            marker = list(color = 'rgb(118, 215, 196)',
            line = list(color = 'rgb(14, 102, 85)',
                        width = 1.5))) %>% 
  add_bars()
fig %>% 
  layout(title = list(text="Puntaje del audio", y = 0.95),
         xaxis = list(title = 'Rating'),
         yaxis = list(title = 'Frecuencia'))

En el siguiente grafico se puede observar como se distribuye el puntaje por especie de ave. Para todas las ves, la media del puntaje es 4, y todas, a excepción de skylar, tiene mas del 50% de los datos con un puntaje igual o superior a 4. Para el caso de la especie skylar, alrededor del 75% de los audios, tienen una puntuación igual o inferior a 4, por lo que está ave podría ser la que tenga desventaja en el modelo, al momento de realizar el entrenamiento, pero esto se verificará al entrenar el modelo y evaluar los resultados.

library(forcats)

fig <- sub_data %>% 
  # mutate(rating = fct_reorder(rating, n, .desc = T)) %>% 
    plot_ly(x=~rating, type = "box", color=~etiqueta) 
fig %>% 
  layout(title = list(text="Puntaje del audio", y = 0.95),
         xaxis = list(title = 'Rating'),
         yaxis = list(title = 'Frecuencia'))

Tiempo

Esta variable representa la hora del día en la que se registró la grabación. Con base en esta variable, se creará una variable categórica que indique el momento del día en que se tomo la grabación: madrugada, mañana, tarde, noche.

library(lubridate)
sub_data$hora <- hour(hm(sub_data$tiempo))
sub_data <- sub_data %>% 
  mutate(etapaDia = case_when(
                              hora >= 0 & hora <=6 ~"Madrugada",
                              hora > 6 & hora <= 12 ~"Mañana",
                              hora > 12 & hora <= 18 ~"Tarde",
                              hora > 18~"Noche",
                            T ~"No especifica"))

fig <- sub_data %>% 
  group_by(etiqueta, etapaDia) %>% 
  summarise(conteo=n()) %>% 
  plot_ly(x = ~etapaDia, y = ~conteo, color = ~etiqueta) %>% 
    layout(title = list(text="Momento del dia en que se tomo el audio", y = 0.95),
         xaxis = list(title = 'Momento del dia'),
         yaxis = list(title = 'Frecuencia'))

fig

Se logra apreciar que para el caso de la especie brnnowl, la mayoría de las grabaciones fueron tomadas en horas de la madrugada, seguido por la noche, lo que es acorde al hecho a que está especie es nocturna. Similar ocurre con la especie comsan, que tambien tiene mayor registro de audio en esos momentos del día, aunque esta especie no es nocturna, por lo que se esperaría mayor registro de audio en momentos como la mañana en lugar de la noche. Para las demás especie, el momento del día en que mayor cantidad de audio se registraron fue en la mañana, seguido de la madrugada.

Distribución Geográfica (MAPA)

Al tener el registro de latitud y longitud de los audios, se puede plasmar un mapa que ayude a observar la distribución geográfica de las especies.

# Mapa de distribución geográfica de las especies

fig <- plot_geo(sub_data, lat = ~latitud, lon = ~longitud)
fig <- fig %>% add_markers(
    text = ~paste(paste("Nombre Común:", Nombre_comun), paste("Nombre Científico:", Nombre_cientifico), sep = "<br />"),
    color = ~etiqueta, symbol = I("circle"), size = I(8), hoverinfo = "text", colors = "viridis"
  ) %>% 
  layout(showlegend=T,
          title = 'Distribución geográfica de las aves', 
         geo = list(showland= T, landcolor = toRGB("gray95"),
                    showcountries = T, subunitcolor = toRGB("gray85"),
                    countrycolor = toRGB("gray85"),
                    countrywidth = 0.5,
                    subunitwidth = 0.5))
fig

• En el mapa se puede observar, que la especie norcar se encuentra distribuida alrededor de México, Estados Unidos y Hawaii.

• La especie skylar, se encuentra distribuida principalmente por Europa, con una presencia reducida en algunos países asiáticos, Australia y en algunos puntos de Estados Unidos

• Para el caso de la especie mallar3, se encuentra distribuida principalmente por Europa, con una presencia reducida en algunos países asiáticos, como Japón, Nueva Zelanda, Estados Unidos y Colombia.

• La especie houspa tiene presencia en varios puntos de todo el continente Americano, también en Europa, parte de Asia, el continente Africano, Australia y Nueva Zelanda.

• La especie comsan se puede observar en gran parte de Europa, algunos puntos de África y Asia.

• Finalmente, la especie brnowl, se encuentra principalmente en Europa, con presencia en algunos puntos del continente Americano, el continente Africano y en algunos países de Asia.

Type (Tipo de sonido)

Esta variable contiene detalles sobre el audio y y el ave presente en el. Se puede encontrar información como el tipo de vocalización emitida por el ave (canto o llamado), el sexo del ave, o en que etapa se encuentra (adulto o juvenil), entre otra información referente a los sonidos que se pueden escuchar en el audio.

# limpieza de la variable type
# se eliminan algunos símbolos, se coloque todo en minúscula y luego se crea una lista con todos los string en la variable
# luego se separan los string para poder obtener los valores único de la columna y poder identificar palabras mas frecuentes
sub_data$type <- str_replace_all(sub_data$type, "[\\[\\]'\\s]", "") %>% 
  str_to_lower()

unique_values <- str_c(sub_data$type, collapse = ",") %>% 
  strsplit(',')
# unique_values[[1]] %>% unique()
head(sub_data) %>% kable()

# creacción de variable sexo
sub_data%>% 
  mutate(sexo=case_when(str_detect(type, "female") & str_detect(type, ",male") ~ "Ambos",
                        str_detect(type, "female")~"Hembra",
                        str_detect(type, ",male")~"Macho",
                        T~"No Especifica")) -> sub_data
# grafico de barras
sub_data %>%
  group_by(sexo) %>% 
  summarise(conteo=n()) %>% 
  plot_ly(x= ~conteo, y=~sexo,  orientation ="h" ) %>% 
  add_bars( marker = list(color = 'rgb(118, 215, 196)',
                      line = list(color = 'rgb(14, 102, 85)',
                                  width = 1.5))) %>% 
      layout(title = list(text="Sexo del Ave"),
         xaxis = list(title = 'Sexo'),
         yaxis = list(title = 'Frecuencia'))

# grafico de barras
sub_data %>%
  group_by(etiqueta, sexo) %>% 
  summarise(conteo=n()) %>% 
  plot_ly(x= ~conteo, y=~sexo,  orientation ="h", color = ~etiqueta) %>% 
  add_bars() %>% 
      layout(title = list(text="Sexo del Ave"),
         xaxis = list(title = 'Sexo'),
         yaxis = list(title = 'Frecuencia'))

sub_data %>% 
  mutate(edad=case_when(str_detect(type, "adult")~"adulto",
                        str_detect(type, "juvenile")~"juvenil",
                        T~"No Especifica")) -> sub_data
sub_data %>%
  group_by(edad) %>% 
  summarise(conteo=n()) %>% 
  plot_ly(x= ~conteo, y=~edad,  orientation ="h" ) %>% 
  add_bars( marker = list(color = 'rgb(118, 215, 196)',
                      line = list(color = 'rgb(14, 102, 85)',
                                  width = 1.5))) %>% 
      layout(title = list(text="Edad del Ave"),
         xaxis = list(title = 'Edad'),
         yaxis = list(title = 'Frecuencia'))

# grafico de barras
sub_data %>%
  group_by(etiqueta, edad) %>% 
  summarise(conteo=n()) %>% 
  plot_ly(x= ~conteo, y=~edad,  orientation ="h", color = ~etiqueta) %>% 
  add_bars() %>% 
      layout(title = list(text="Edad del Ave"),
         xaxis = list(title = 'Edad'),
         yaxis = list(title = 'Frecuencia'))

sub_data %>% 
  mutate(tipo=case_when(str_detect(type, "call")~"llamado",
                        str_detect(type, "song")~"canto",
                        T~"No Especifica")) -> sub_data
sub_data %>%
  group_by(tipo) %>% 
  summarise(conteo=n()) %>% 
  plot_ly(x= ~conteo, y=~tipo,  orientation ="h" ) %>% 
  add_bars( marker = list(color = 'rgb(118, 215, 196)',
                      line = list(color = 'rgb(14, 102, 85)',
                                  width = 1.5))) %>% 
      layout(title = list(text="Tipo de vocalización del Ave"),
         xaxis = list(title = 'Vocalización'),
         yaxis = list(title = 'Frecuencia'))

# grafico de barras
sub_data %>%
  group_by(etiqueta, tipo) %>% 
  summarise(conteo=n()) %>% 
  plot_ly(x= ~conteo, y=~tipo,  orientation ="h", color = ~etiqueta) %>% 
  add_bars() %>% 
      layout(title = list(text="Tipo de vocalización del Ave"),
         xaxis = list(title = 'Tipo'),
         yaxis = list(title = 'Frecuencia'))

CNN

El primer modelo que se entrenará es una red neuronal convolucional (CNN). La arquitectura propuesta para este red neuronal es la siguiente:

# Definir el modelo
model <- keras_model_sequential()
weight_decay = 1e-4

# Añadir capas al modelo
model %>%
  
  layer_conv_2d(filters = 32, kernel_size = 3, strides = 1, padding = 'same',
                activation = 'relu', kernel_regularizer = regularizer_l2(weight_decay), 
                input_shape = c(target_size, 3) ) %>%
  
  layer_max_pooling_2d(pool_size = c(2, 2)) %>%
  
  layer_dropout(rate = 0.2) %>%
  
  layer_conv_2d(filters = 64, kernel_size = 3, padding = 'same', activation = 'relu',
                kernel_regularizer = regularizer_l2(weight_decay)) %>%
  
  layer_max_pooling_2d(pool_size = c(2, 2)) %>%
  
  layer_dropout(rate = 0.3) %>%
  
  layer_conv_2d(filters = 128, kernel_size = 3, padding = 'same', activation = 'relu',
                kernel_regularizer = regularizer_l2(weight_decay)) %>%
  
  layer_max_pooling_2d(pool_size = c(2, 2)) %>%
  
  layer_dropout(rate = 0.4) %>%
  
  layer_flatten() %>%
  layer_dense(units = 256, activation = 'relu') %>%
  layer_dense(units = 6, activation = 'softmax')

model

## Model: "sequential"
## ________________________________________________________________________________
##  Layer (type)                       Output Shape                    Param #     
## ================================================================================
##  conv2d_2 (Conv2D)                  (None, 224, 224, 32)            896         
##  max_pooling2d_2 (MaxPooling2D)     (None, 112, 112, 32)            0           
##  dropout_2 (Dropout)                (None, 112, 112, 32)            0           
##  conv2d_1 (Conv2D)                  (None, 112, 112, 64)            18496       
##  max_pooling2d_1 (MaxPooling2D)     (None, 56, 56, 64)              0           
##  dropout_1 (Dropout)                (None, 56, 56, 64)              0           
##  conv2d (Conv2D)                    (None, 56, 56, 128)             73856       
##  max_pooling2d (MaxPooling2D)       (None, 28, 28, 128)             0           
##  dropout (Dropout)                  (None, 28, 28, 128)             0           
##  flatten (Flatten)                  (None, 100352)                  0           
##  dense_1 (Dense)                    (None, 256)                     25690368    
##  dense (Dense)                      (None, 6)                       1542        
## ================================================================================
## Total params: 25,785,158
## Trainable params: 25,785,158
## Non-trainable params: 0
## ________________________________________________________________________________

Para entrenar el modelo se considera la función de perdida de categorical crossentropy, y una tasa de aprendizaje de 0.01.

model %>% 
  compile(
    loss = "categorical_crossentropy",
    optimizer = optimizer_adam(lr = 0.01),
    metrics = c("accuracy", metric_auc())
  )

# Ajuste del modelo
history <- model %>% 
  fit(
  # training data
  train_image_array_gen,

  # training epochs
  steps_per_epoch = as.integer(train_samples / batch_size), 
  epochs = 10, 
  
  # validation data
  validation_data = val_image_array_gen,
  validation_steps = as.integer(valid_samples / batch_size)
)

# save_model_tf(model, "cnn_aves_1")
# write.csv(history, "historyCNN_1.csv")

Se carga el modelo previamente entrenado para evitar variaciones en los resultados:

model <- load_model_tf("C:\\Users\\kaes1\\Documents\\DataViz\\App_Aves\\cnn_aves_1")
historyCNN1 <- read.csv("C:\\Users\\kaes1\\Documents\\DataViz\\App_Aves\\APP_AVES\\historyCNN_1.csv")

Del proceso de entrenamiento se extrae la métrica del AUC para cada una de las épocas de entrenamiento, y para los datos del train y de validación. Tambien se tiene información de la perdida del modelo, y de la métrica de exactitud o Accuracy.

 # grafica de perdida train vs val
loss_df <- historyCNN1[historyCNN1$metric=="loss", ]

fig1 <- loss_df %>% 
 plot_ly(x=~epoch, y=~value, color = ~data, type="scatter", mode = 'lines+markers') %>% 
 layout(legend = list(orientation = 'v', y = 1, x=0.5), 
        title = "Loss - Train vs Val")
fig1

 # grafica de auc train vs val

auc_df <- historyCNN1[historyCNN1$metric=="auc_2", ]

fig2 <- auc_df %>% 
 plot_ly(x=~epoch, y=~value, color = ~data, type="scatter", mode = 'lines+markers') %>% 
 layout(legend = list(orientation = 'v', y = 0.1, x=0.5), 
        title = "AUC - Train vs Val")
fig2   

 # Graficar Accuracy
acc_df <- historyCNN1[historyCNN1$metric=="accuracy", ]
   
fig3<- acc_df %>% 
 plot_ly(x=~epoch, y=~value, color = ~data, type="scatter", mode = 'lines+markers') %>% 
 layout(legend = list(orientation = 'v', y = 0.1, x=0.5),
        title="Accuracy - Train vs Val")
fig3   

Para hacer el test del modelo, se deben transformar tambien las imágenes de prueba, aplicando el re-dimensionamiento y re-escalando los pixeles de las imágenes.

# funcion para convertir las imagenes en arreglos vectoriales]
image_prep <- function(x) {
  arrays <- lapply(x, function(path) {
    img <- image_load(path, target_size = target_size, 
                      grayscale = F 
                      )
    
    x <- image_to_array(img)
    x <- array_reshape(x, c(1, dim(x)))
    x <- x/255 # rescale image pixel
  })
  do.call(abind::abind, c(arrays, list(along = 1)))
}

test_path <- 'C:/Users/kaes1/Desktop/MachineLearningUN/Proyecto/test/'
test_files <- list()

for (folder in folder_list) {
 test_files <- append(test_files, paste0(paste0(test_path, folder),"/",list.files(paste0(test_path, folder))))
}

test_data <- data.frame(file_name = test_files %>% as.character()) %>% 
  mutate(class = str_extract(test_files %>% as.character(), "brnowl|comsan|houspa|mallar3|norcar|skylar"))

head(test_data, 5) %>% kable()

test_x <- image_prep(test_data$file_name)
dim(test_x)

## [1] 600 224 224   3

mtx_test <- model %>% predict(test_x) 

## 19/19 - 8s - 8s/epoch - 395ms/step

head(mtx_test, 10) %>% kable()

pred_test <- mtx_test %>% k_argmax()
prob_test <- apply(mtx_test, 1, max, na.rm=TRUE)
head(pred_test, 10)

## tf.Tensor([4 0 3 0 0 4 0 0 4 0], shape=(10), dtype=int64)

# Convertir encoding a etiqueta
decode <- function(x){
  case_when(x == 0 ~ "brnowl",
            x == 1 ~ "comsan",
            x == 2 ~ "houspa",
            x == 3 ~ "mallar3",
            x == 4 ~ "norcar",
            x == 5 ~ "skylar"
            
            )
}

pred_test <- sapply(pred_test, decode) 

head(pred_test, 10)

##  [1] "norcar"  "brnowl"  "mallar3" "brnowl"  "brnowl"  "norcar"  "brnowl" 
##  [8] "brnowl"  "norcar"  "brnowl"

confusionMatrix(as.factor(pred_test), 
                as.factor(test_data$class)
                ) -> cmatrix
cmatrix

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction brnowl comsan houspa mallar3 norcar skylar
##    brnowl      69      0      3      12      3      0
##    comsan       7     91      3       8      5     18
##    houspa       6      4     86       2     10     24
##    mallar3      7      1      0      65      0      1
##    norcar      10      2      6      13     82      6
##    skylar       1      2      2       0      0     51
## 
## Overall Statistics
##                                                
##                Accuracy : 0.74                 
##                  95% CI : (0.7029, 0.7747)     
##     No Information Rate : 0.1667               
##     P-Value [Acc > NIR] : < 0.00000000000000022
##                                                
##                   Kappa : 0.688                
##                                                
##  Mcnemar's Test P-Value : 0.0000000004847      
## 
## Statistics by Class:
## 
##                      Class: brnowl Class: comsan Class: houspa Class: mallar3
## Sensitivity                 0.6900        0.9100        0.8600         0.6500
## Specificity                 0.9640        0.9180        0.9080         0.9820
## Pos Pred Value              0.7931        0.6894        0.6515         0.8784
## Neg Pred Value              0.9396        0.9808        0.9701         0.9335
## Prevalence                  0.1667        0.1667        0.1667         0.1667
## Detection Rate              0.1150        0.1517        0.1433         0.1083
## Detection Prevalence        0.1450        0.2200        0.2200         0.1233
## Balanced Accuracy           0.8270        0.9140        0.8840         0.8160
##                      Class: norcar Class: skylar
## Sensitivity                 0.8200       0.51000
## Specificity                 0.9260       0.99000
## Pos Pred Value              0.6891       0.91071
## Neg Pred Value              0.9626       0.90993
## Prevalence                  0.1667       0.16667
## Detection Rate              0.1367       0.08500
## Detection Prevalence        0.1983       0.09333
## Balanced Accuracy           0.8730       0.75000

library(tibble)
library(tidyr)
cmatrix$byClass %>% as.data.frame() -> metrics_pcf
metrics_pcf %>% rownames_to_column() -> metrics_pcf
metrics_pcf$rowname <- str_replace(metrics_pcf$rowname, "Class: ", "")
metrics_pcf <- metrics_pcf %>% select(rowname, Precision, Recall, F1)

metrics_pcf <- gather(metrics_pcf,
                 "Precision", "Recall", "F1",
                 key = "Metrica",
                 value = "Score")
metrics_pcf

##    rowname   Metrica     Score
## 1   brnowl Precision 0.7931034
## 2   comsan Precision 0.6893939
## 3   houspa Precision 0.6515152
## 4  mallar3 Precision 0.8783784
## 5   norcar Precision 0.6890756
## 6   skylar Precision 0.9107143
## 7   brnowl    Recall 0.6900000
## 8   comsan    Recall 0.9100000
## 9   houspa    Recall 0.8600000
## 10 mallar3    Recall 0.6500000
## 11  norcar    Recall 0.8200000
## 12  skylar    Recall 0.5100000
## 13  brnowl        F1 0.7379679
## 14  comsan        F1 0.7844828
## 15  houspa        F1 0.7413793
## 16 mallar3        F1 0.7471264
## 17  norcar        F1 0.7488584
## 18  skylar        F1 0.6538462

fig <- plot_ly(metrics_pcf, x = ~Metrica, y = ~Score, color = ~rowname) 
fig %>% add_markers(size=20) %>% layout(title="Metricas de desempeño")

cmatrix$table %>% as.data.frame() %>%  plot_ly(x = ~Prediction, y =~Reference, z= ~Freq, type = "heatmap") %>% 
  layout(title="Matriz de Confusión")-> fig
fig

library(fastDummies)
library(yardstick)
mtx_test %>% as.data.frame() -> y_scores
names(y_scores) <- paste0("pred_", birds)

y_onehot <- dummy_cols(test_data$class)
colnames(y_onehot) <- c("drop", birds)
y_onehot <- subset(y_onehot, select = -c(drop))

z <- cbind(y_scores, y_onehot)

z$brnowl <- as.factor(z$brnowl)
roc_brnowl <- roc_curve(data = z, brnowl, pred_brnowl)
roc_brnowl$specificity <- 1 - roc_brnowl$specificity
colnames(roc_brnowl) <- c('threshold', 'tpr', 'fpr')
auc_brnowl <- roc_auc(data = z, brnowl, pred_brnowl)
auc_brnowl <- auc_brnowl$.estimate
brnowl <- paste('brnowl (AUC=',toString(round(1-auc_brnowl,2)),')',sep = '')


z$houspa <- as.factor(z$houspa)
roc_houspa <- roc_curve(data = z, houspa, pred_houspa)
roc_houspa$specificity <- 1 - roc_houspa$specificity
colnames(roc_houspa) <- c('threshold', 'tpr', 'fpr')
auc_houspa <- roc_auc(data = z, houspa, pred_houspa)
auc_houspa <- auc_houspa$.estimate
houspa <- paste('houspa (AUC=',toString(round(1-auc_houspa,2)),')',sep = '')

z$comsan <- as.factor(z$comsan)
roc_comsan <- roc_curve(data = z, comsan, pred_comsan)
roc_comsan$specificity <- 1 - roc_comsan$specificity
colnames(roc_comsan) <- c('threshold', 'tpr', 'fpr')
auc_comsan <- roc_auc(data = z, comsan, pred_comsan)
auc_comsan <- auc_comsan$.estimate
comsan <- paste('comsan (AUC=',toString(round(1-auc_comsan,2)),')',sep = '')

z$mallar3 <- as.factor(z$mallar3)
roc_mallar3 <- roc_curve(data = z, mallar3, pred_mallar3)
roc_mallar3$specificity <- 1 - roc_mallar3$specificity
colnames(roc_mallar3) <- c('threshold', 'tpr', 'fpr')
auc_mallar3 <- roc_auc(data = z, mallar3, pred_mallar3)
auc_mallar3 <- auc_mallar3$.estimate
mallar3 <- paste('mallar3 (AUC=',toString(round(1-auc_mallar3,2)),')',sep = '')

z$norcar <- as.factor(z$norcar)
roc_norcar <- roc_curve(data = z, norcar, pred_norcar)
roc_norcar$specificity <- 1 - roc_norcar$specificity
colnames(roc_norcar) <- c('threshold', 'tpr', 'fpr')
auc_norcar <- roc_auc(data = z, norcar, pred_norcar)
auc_norcar <- auc_norcar$.estimate
norcar <- paste('norcar (AUC=',toString(round(1-auc_norcar,2)),')',sep = '')

z$skylar <- as.factor(z$skylar)
roc_skylar <- roc_curve(data = z, skylar, pred_skylar)
roc_skylar$specificity <- 1 - roc_skylar$specificity
colnames(roc_skylar) <- c('threshold', 'tpr', 'fpr')
auc_skylar <- roc_auc(data = z, skylar, pred_skylar)
auc_skylar <- auc_skylar$.estimate
skylar <- paste('skylar (AUC=',toString(round(1-auc_skylar,2)),')',sep = '')

fig <- plot_ly()%>%
  add_segments(x = 0, xend = 1, y = 0, yend = 1, line = list(dash = "dash", color = 'black'), showlegend = FALSE) %>%
  add_trace(data = roc_brnowl,x = ~fpr, y = ~tpr, mode = 'lines', name = brnowl, type = 'scatter')%>%
  add_trace(data = roc_comsan,x = ~fpr, y = ~tpr, mode = 'lines', name = comsan, type = 'scatter')%>%
  add_trace(data = roc_houspa,x = ~fpr, y = ~tpr, mode = 'lines', name = houspa, type = 'scatter')%>%
    add_trace(data = roc_mallar3,x = ~fpr, y = ~tpr, mode = 'lines', name = mallar3, type = 'scatter')%>%
    add_trace(data = roc_norcar,x = ~fpr, y = ~tpr, mode = 'lines', name = norcar, type = 'scatter')%>%
    add_trace(data = roc_skylar,x = ~fpr, y = ~tpr, mode = 'lines', name =skylar, type = 'scatter')%>%
  layout(xaxis = list(
    title = "False Positive Rate"
  ), yaxis = list(
    title = "True Positive Rate"
  ),legend = list(x = 0.8, y = 0.1))
fig

CNN-SVC

El segundo modelo a evaluar, es un híbrido entre CNN y Support Vector Classifier. Este modelo utilizará la red convolucional entrenada anteriormente para extraer atributos de las imágenes, teniendo en cuenta la penúltima capa de la red.

Se crea un data frame con las rutas de las imágenes y la etiqueta correspondiente a la especie para cada una.

train_path <- 'C:/Users/kaes1/Desktop/MachineLearningUN/Proyecto/train/'
train_files <- list()

for (folder in folder_list) {
 train_files <- append(train_files, paste0(paste0(train_path, folder),"/",list.files(paste0(train_path, folder))))
}

train_data <- data.frame(file_name = train_files %>% as.character()) %>% 
  mutate(class = str_extract(train_files %>% as.character(), "brnowl|comsan|houspa|mallar3|norcar|skylar"))

head(train_data, 5) %>% kable()

Debido a limitantes en cuanto a la memoria del equipo de computo utilizado, para entrenar este modelo híbrido, se considera un tamaño de entrenamiento menor al utilizado en el modelo de CNN. En este caso el set de entrenamiento tiene un tamaño de 1200 imágenes.

train_data_ <-  createDataPartition(y =train_data$class, p = 0.5, times = 1, list =F )
train_data_<- train_data[train_data_,]
dim(train_data_)

## [1] 1200    2

train_x <- image_prep(train_data_$file_name)
dim(train_x)

## [1] 1200  224  224    3

Se crea una función que permite hacer uso del modelo CNN, y de la capa número 11 (la capa densa previa a la clasificación) para extraer los atributos de las imágenes de entrenamiento.

feature_extractor <- keras_model(
  inputs = model$inputs,
  outputs = get_layer(model, index = 11)$output
)

# extraer atributos de los datos de entrenamiento
features <- feature_extractor(train_x)

# features[[11]] %>% as.matrix() %>% as.data.frame() -> data_prueba
features %>% as.matrix() %>% as.data.frame() -> data_train
data_train <- data_train %>% mutate("Etiqueta" = train_data_$class %>% as.factor())
data_train %>% head() %>%  kable()

En la capa número 11 de la CNN, se tiene como salida 256 características extraídas de las imágenes, las cuales serán nuestros datos de entrenamiento para el modelo de SVC (SVM).

dim(data_train)

## [1] 1200  257

Para los datos de prueba se aplica el mismo proceso, con el fin de que se tenga la misma estructura que los datos de entrenamiento.

features_test <- feature_extractor(test_x)
features_test %>% as.matrix() %>% as.data.frame() -> data_test
data_test <- data_test %>% mutate("Etiqueta" = test_data$class %>% as.factor())
data_test %>% head() %>% kable()