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'))

Entrenamiento

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)))
}

Predicción con el test set

Se obtiene el listado de imágenes del set de prueba:

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()

Ahora se procesan las imágenes para que tengan la estructura requerida por el modelo de CNN. Se verifica que los datos si tengan la estructura especificada.

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

## [1] 600 224 224   3

Se procede a hacer la predicción con los datos de prueba, y se obtiene una matriz con las probabilidades de clasificación para cada una de las especies.

mtx_test <- model %>% predict(test_x) 

## 19/19 - 9s - 9s/epoch - 498ms/step

head(mtx_test, 5) 

##            [,1]            [,2]          [,3]       [,4]      [,5]
## [1,] 0.21947034 0.0000008590754 0.00003032223 0.02016852 0.7603290
## [2,] 0.48916000 0.0072226556949 0.06774273515 0.22360452 0.2108659
## [3,] 0.08020142 0.1564177870750 0.08535960317 0.31130502 0.2898542
## [4,] 0.78099644 0.0099687641487 0.00296757231 0.01959453 0.1830159
## [5,] 0.55572104 0.0347285307944 0.00390776759 0.21312954 0.1858690
##                [,6]
## [1,] 0.000000960562
## [2,] 0.001404142473
## [3,] 0.076861850917
## [4,] 0.003456914797
## [5,] 0.006644142326

Ahora se obtiene la clasificación final para cada registro, obteniendo la posición en la matriz de la probabilidad más grande.

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)

Para facilitar la compresión de los resultados, se hace la decodificación de la variable respuesta, para poder identificar a que ave se hace referencia.

# 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"

Matrix de Confusión

A continuación se obtiene la matriz de confusión para predicción realizada. En los resultados obtenidos, se observa la matriz de confusión, junto con otras métricas de desempeño, como la especificidad, sensibilidad, Accuracy, etc.

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

En el siguiente gráfico se observan las métricas de F1-Score, Precision y Recall. Se logra ver que los valores de Recall para cada especie son los que tienen mayor dispersión, siendo la especie skylar la que tiene el valor mas bajo, lo que indica que el modelo solo logró identificar correctamente el 51% de los espectrogramas pertenecientes a esta especie. Por otro lado, esta especie tiene una Precision alta, ya que el 91% de los espectrogramas clasificados en esa especie, si pertenencia a esa especie. Similar ocurre con las especies brnowl y mallar3

Lo contrario sucede con la especie houspa, la cual tiene un Recall, 86%, pero una Precision mas baja, 65%. Similar a esto ocurre con la especie comsan y norcar.

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")

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

Se puede ver en la matriz, que las especies con menor número de clasificaciones correctas son brnowl, mallar3 y skylar. La especie con mayor número de clasificaciones correctas es comsan, seguida de houspa y norcar. Esto podría indicar que el modelo tiene un sesgo hacia la especie estas tres especies, y no es tan bueno para identificar las demás especies.

Curva ROC

A continuación se traza la curva ROC para cada una de las especies, bajo la modalidad de One vs All, donde se mira que tan bueno es el modelo para identificar cada especie individualmente, frente a las demás.

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

Se obtienen las etiqueta de respuesta con one hot encoding.

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

Ahora se gráfica para cada especie la Curva ROC, junto con el área bajo la curva (AUC), correspondiente.

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 = '')

En general, los valores de AUC se encuentran todos por encima de 0.7, lo que indica que el modelo es bueno para distinguir a las especies individualmente. Sin embargo, los análisis realizados con las otras métricas indican que el modelo podría tener un sesgo hacia algunas especies, lo cual se debería confirmar en analisis mas profundos de los resultados del modelo.

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