Análisis de una base de datos de biodiversidad: avifauna en la Reserva Thomas vander Hammen

Introducción

El presente notebook muestra el código utilizado para realizar el estudio presentado en el artículo “Metodologías de estudio estadístico de bases de datos de biodiversidad: el caso de una base de datos de la biodiversidad de avifauna en la reserva Thomas van der Hammen.” y tiene como objetivo detallar los procedimientos y análisis estadísticos aplicados a la base de datos estudiada. Se describen los pasos seguidos desde la limpieza y preparación de los datos hasta la aplicación de técnicas estadísticas avanzadas y la visualización de resultados.

Preparación de la base de datos

Es importante destacar que la base de datos original (Barrera, 2023) contenía 53 variables. Antes de trabajar en el entorno de R, se realizó un primer filtro de información, reduciendo el número de variables a 18, las cuales formaron parte del estudio final. Esta versión de la base de datos puede descargarse ingresando a este enlace.

Cargue de la base de datos

library(readr)
dbAves=read_delim("https://drive.google.com/uc?id=1SvGS1WHzoLafehGMnvDPGxGvyQLbgRKd",
                     delim = ",",
                     escape_double = FALSE,
                     trim_ws = TRUE)

Rows: 3499 Columns: 18
── Column specification ───────────────────────────────────────────
Delimiter: ","
chr (16): id, sex, behavior, associatedTaxa, occurrenceRemarks,...
dbl  (2): decimalLatitude, decimalLongitude

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

summary(dbAves)

      id                sex              behavior        
 Length:3499        Length:3499        Length:3499       
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                         
                                                         
                                                         
 associatedTaxa     occurrenceRemarks    habitat         
 Length:3499        Length:3499        Length:3499       
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                         
                                                         
                                                         
 verbatimLocality   decimalLatitude   decimalLongitude   
 Length:3499        Min.   : 476978   Min.   :-74101169  
 Class :character   1st Qu.:4767098   1st Qu.:-74100289  
 Mode  :character   Median :4767918   Median :-74098822  
                    Mean   :4551481   Mean   :-72875474  
                    3rd Qu.:4769195   3rd Qu.:-74097725  
                    Max.   :4772912   Max.   : -7410075  
 verbatimLatitude   verbatimLongitude  scientificName    
 Length:3499        Length:3499        Length:3499       
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                         
                                                         
                                                         
    order              family             genus          
 Length:3499        Length:3499        Length:3499       
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                         
                                                         
                                                         
 parentEventID      estructuraSocial   alturaObservacion 
 Length:3499        Length:3499        Length:3499       
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                         
                                                         
                                                         

Primera fase

Del resumen obtenido anteriormente se observa que:

• Aunque las variables decimalLatitude y decimalLongitude se identifican como cuantitativas, no son correctas, ya que sus partes decimales se convirtieron en enteros. Para corregir esto, se utilizarán las variables verbatimLatitude y verbatimLongitude. Una vez realizada esta modificación, estas dos últimas variables pueden eliminarse, ya que se habrá obtenido la medida numérica correcta.
• Las variables que R clasifica como character deben convertirse a factores.
• La variable alturaObservación debe ser tratada para que tome valores numéricos.

Se inicia con los procedimientos para abordar estos hallazgos.

Variables decimalLatitude y decimalLongitude

Para hacer los cambios en estas variables, se utilizarán las variables verbatimLatitude y verbatimLongitude para realizar la conversión correcta de las coordenadas de su forma GMS a decimal. Se programará una función que tome la cadena de texto con las coordenadas GMS, identifique sus componentes, los convierta a valores numéricos y luego realice las operaciones aritméticas necesarias para obtener el valor decimal de la coordenada. Una vez realizada la conversión, se eliminarán las variables verbatimLatitude y verbatimLongitude.

convLat <- function(ejemplo) {
  partes1 =strsplit(ejemplo, "\u00B0")
  grados = as.numeric(partes1[[1]][1])
  partes2=strsplit(partes1[[1]][2], "\u00B4")
  minutos =as.numeric(partes2[[1]][1])
  segundos =as.numeric(partes2[[1]][2])
  decimalLat = grados + minutos / 60 + segundos / 3600
  return(decimalLat)}
convLon=function(ejemplo){
  partes1 =strsplit(ejemplo, "\u00B0")
  grados = as.numeric(partes1[[1]][1])
  partes2=strsplit(partes1[[1]][2], "\u00B4")
  minutos =as.numeric(partes2[[1]][1])
  segundos =as.numeric(partes2[[1]][2])
  decimalLon = (grados + minutos / 60 + segundos / 3600)*-1
  return(decimalLon)
}
dbAves$decimalLatitude=sapply(dbAves$verbatimLatitude,convLat)
dbAves$decimalLongitude=sapply(dbAves$verbatimLongitude,convLon)
dbAves=subset(dbAves,select = -verbatimLatitude)
dbAves=subset(dbAves,select = -verbatimLongitude)

Conversión de variables tipo character o carácter.

Para hacer un correcto tratamiento de variables cualitativas en el entorno de R, es necesario que estas sean de tipo factor. En este caso, se realizará la conversión de las variables que inicialmente se identificaron como de tipo carácter (incluida la variable alturaObservacion para poder explorarla, a pesar de que debería ser una variable numérica).

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

dbAves=dbAves %>%
  mutate(across(where(is.character), as.factor))

Variable alturaObservacion

En este punto, esta variable es de tipo factor, por lo que un buen primer paso es revisar las categorías que la componen.

unique(dbAves$alturaObservacion)

 [1] #N/D  10    5     3     4     2     6     7     9     8    
[11] <0.6  15    12    13    1.5-2 1     1.5   11    0.6   4.5  
[21] 2.5   Bajo  20    >4    2.4   0.80  0.8   3-5  
28 Levels: #N/D <0.6 >4 0.6 0.8 0.80 1 1.5 1.5-2 10 11 12 ... Bajo

Como se puede evidenciar, hay categorías que no son propiamente un número, como “<0.6”, “1.5-2”, “Bajo”, “>4” y “3-5”. Con respecto a la categoría “#N/D”, son datos de los cuales no se tiene un valor registrado, por lo que no se deben modificar y deben tratarse como datos faltantes.

Dado que esta es una variable de naturaleza cuantitativa, sería interesante poder tratarla como tal. Por lo tanto, para este estudio se realizarán modificaciones a estas categorías y a los registros que las contengan. Estas modificaciones se realizarán teniendo en cuenta los siguientes parámetros:

• La categoría “#N/D” debe modificarse para que tome valores de NA.
• A las categorías que denoten intervalos se les asignará el valor medio entre los límites del mismo.
• A la categoría “<0.6” se le designará el valor de 0.5.
• A la categoría “>4” se le asignará el valor de 5.
• A la categoría “Bajo” se le designará el valor de 0.3.

dbAves <- dbAves %>%
  mutate(alturaObservacion = case_when(
    alturaObservacion == "#N/D" ~ NA_real_,
    alturaObservacion == "1.5-2" ~ 1.7,
    alturaObservacion == "3-5" ~ 4,
    alturaObservacion == "<0.6" ~ 0.5,
    alturaObservacion == ">4" ~ 5,
    alturaObservacion == "Bajo" ~ 0.3,
    TRUE ~ as.numeric(as.character(alturaObservacion))
  ))

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `alturaObservacion = case_when(...)`.
Caused by warning:
! NAs introducidos por coerción

Para finalizar esta primera fase de limpieza de la base de datos, la variable ID se convierte en el identificador de fila y no se tomará en cuenta como una variable a partir de este momento.

rownames(dbAves)=dbAves$id

Warning: Setting row names on a tibble is deprecated.

dbAves$id=NULL

summary(dbAves)

          sex                                       behavior   
 Desconocido:1666   Vuelo                               :1860  
 Hembra     :  66   Descanso                            : 842  
 Juvenil    :   5   Movimiento entre vegetación         : 385  
 Macho      : 123   Alimentación                        : 130  
 NA's       :1639   Alimentación | Movimiento vegetación:  77  
                    Autoacicalamiento                   :  53  
                    (Other)                             : 152  
                          associatedTaxa
 Huésped: Cenchrus clandestinus  : 236  
 Huésped: Smallanthus pyramidalis: 202  
 Huésped: Alnus acuminata        :  82  
 Huésped: Fraxinus uhdei         :  61  
 Huésped: Bacharis latifolia     :  52  
 (Other)                         : 491  
 NA's                            :2375  
           occurrenceRemarks                         habitat    
 Monitoreo 5, mañana: 520    Bosque mixto nativo- exótica: 365  
 Monitoreo 5, tarde : 430    Bosque Nativo               : 710  
 Monitoreo 4, mañana: 403    Espejo de agua              :  14  
 Monitoreo 3, tarde : 378    Matorral pasto Kikuyo       : 288  
 Monitoreo 3, mañana: 348    NA's                        :2122  
 Monitoreo 1, tarde : 285                                       
 (Other)            :1135                                       
              verbatimLocality decimalLatitude decimalLongitude
 Hacienda la Conejera : 467    Min.   :4.765   Min.   :-74.10  
 Hacienda las Mercedes:3032    1st Qu.:4.767   1st Qu.:-74.10  
                               Median :4.768   Median :-74.10  
                               Mean   :4.769   Mean   :-74.10  
                               3rd Qu.:4.770   3rd Qu.:-74.10  
                               Max.   :4.773   Max.   :-74.07  
                                                               
               scientificName             order     
 Zenaida auriculata   : 611   Passeriformes  :2407  
 Zonotrichia capensis : 492   Columbiformes  : 682  
 Turdus fuscater      : 295   Pelecaniformes : 112  
 Notiochelidon murina : 294   Apodiformes    :  94  
 Molothrus bonariensis: 209   Psittaciformes :  68  
 Spinus psaltria      : 159   Accipitriformes:  61  
 (Other)              :1439   (Other)        :  75  
          family               genus      parentEventID
 Columbidae  : 682   Zenaida      : 611   E:1: 467     
 Emberizidae : 525   Zonotrichia  : 492   E:2: 345     
 Hirundinidae: 419   Notiochelidon: 296   E:3:1280     
 Turdidae    : 303   Turdus       : 295   E:4:1407     
 Icteridae   : 286   Spinus       : 222                
 Thraupidae  : 280   Molothrus    : 209                
 (Other)     :1004   (Other)      :1374                
           estructuraSocial alturaObservacion
 Bandada           : 264    Min.   : 0.300   
 Grupo coespecífico:1222    1st Qu.: 2.000   
 Grupo mixto       :  43    Median : 4.000   
 Pareja            : 809    Mean   : 4.846   
 Solitario         :1161    3rd Qu.: 7.000   
                            Max.   :20.000   
                            NA's   :2532     

head(dbAves)

Segunda fase

Como se puede evidenciar en el último resumen emitido, se observa que hay variables que tienen datos faltantes (aparecen bajo la categoría “NA’s”). Se procede a revisar el porcentaje de datos faltates por variable.

falt=colSums(is.na(dbAves))
porcFalt=round((falt/nrow(dbAves))*100,0)
dfFaltantes=data.frame(Faltantes=falt,
                       Porcentaje=porcFalt)
 dfFaltantes[dfFaltantes$Porcentaje != 0, ]

La presencia de datos faltantes es común al analizar bases de datos de biodiversidad, lo cual puede afectar la confiabilidad de los análisis. En este estudio, se ha decidido eliminar del conjunto de datos aquellas variables que presenten un porcentaje de datos faltantes superior al 65%. En este caso, las variables alturaObservación y associatedTaxa cumplen con este criterio y serán excluidas del análisis.

dbAves=subset(dbAves, select=-alturaObservacion)
dbAves=subset(dbAves, select=-associatedTaxa)
sum(complete.cases(dbAves))

[1] 736

Además, se ha observado que la base de datos contiene un total de 736 registros completos. Por lo tanto, se ha decidido trabajar exclusivamente con estos registros. Se procederá a generar una copia de la base de datos que contenga únicamente estos registros, con el fin de realizar un análisis más detallado de las variables y de la base de datos en general.

sinFalt=dbAves[complete.cases(dbAves),]
write.csv(sinFalt, file = "avesSinFal.csv")

Puede descargarse una copia de esta versión de la base de datos ingresando al siguiente enlace.

Tercera fase

Para esta fase de preparación de los datos, se inicia con la revisión de un resumen de los mismos.

summary(sinFalt)

          sex                                      behavior  
 Desconocido:644   Descanso                            :444  
 Hembra     : 34   Movimiento entre vegetación         :102  
 Juvenil    :  3   Alimentación                        : 52  
 Macho      : 55   Alimentación | Movimiento vegetación: 51  
                   Vuelo                               : 41  
                   Vuelo | Movimiento vegetación       : 11  
                   (Other)                             : 35  
           occurrenceRemarks                         habitat   
 Monitoreo 5, mañana:211     Bosque mixto nativo- exótica:197  
 Monitoreo 5, tarde :139     Bosque Nativo               :442  
 Monitoreo 1, tarde : 74     Espejo de agua              :  3  
 Monitoreo 3, tarde : 44     Matorral pasto Kikuyo       : 94  
 Monitoreo 1, mañana: 43                                       
 Monitoreo 4, mañana: 43                                       
 (Other)            :182                                       
              verbatimLocality decimalLatitude decimalLongitude
 Hacienda la Conejera :175     Min.   :4.765   Min.   :-74.10  
 Hacienda las Mercedes:561     1st Qu.:4.768   1st Qu.:-74.10  
                               Median :4.770   Median :-74.10  
                               Mean   :4.769   Mean   :-74.09  
                               3rd Qu.:4.772   3rd Qu.:-74.10  
                               Max.   :4.773   Max.   :-74.07  
                                                               
               scientificName            order    
 Zonotrichia capensis :137    Passeriformes :571  
 Turdus fuscater      :119    Columbiformes :106  
 Zenaida auriculata   : 66    Apodiformes   : 21  
 Molothrus bonariensis: 40    Psittaciformes: 21  
 Patagioenas fasciata : 40    Pelecaniformes:  6  
 Diglossa humeralis   : 34    Cuculiformes  :  4  
 (Other)              :300    (Other)       :  7  
         family            genus     parentEventID
 Emberizidae:144   Zonotrichia:137   E:1:175      
 Turdidae   :122   Turdus     :119   E:2:190      
 Columbidae :106   Zenaida    : 66   E:3:189      
 Thraupidae : 76   Diglossa   : 44   E:4:182      
 Icteridae  : 70   Molothrus  : 40                
 Tyrannidae : 51   Patagioenas: 40                
 (Other)    :167   (Other)    :290                
           estructuraSocial
 Bandada           :  0    
 Grupo coespecífico:228    
 Grupo mixto       :  7    
 Pareja            :220    
 Solitario         :281    
                           
                           

Como se ha podido observar hasta el momento, el tratamiento de variables cualitativas a menudo requiere un manejo especial de cadenas de texto. A través del resumen anterior, se identificaron otras variables cualitativas en esta base de datos que, debido a la cantidad de categorías o a los valores que algunos registros pueden tomar, podrían estar introduciendo más variabilidad en el conjunto de datos. Estas variables son sex, behavior y occurrenceRemarks. A continuación se presentan los hallazgos y los tratamientos correspondientes.

Variable sex

sinFalt$sex=droplevels(sinFalt$sex)
tabSex=as.data.frame(table(sinFalt$sex))
colnames(tabSex)=c("Sexo","Frecuencia")
tabSex

En el caso de esta variable, aunque hay pocas categorías, la categoría “Juvenil” cuenta con solo tres observaciones. Para reducir el número de categorías de la variable y disminuir su variabilidad, en este estudio se ha decidido clasificar estas tres observaciones como “Desconocido”. Esta elección se basa en que determinar el sexo de un ave a simple vista es factible solo en la mitad de las especies y, usualmente, cuando estas están completamente desarrolladas (es decir, en su etapa adulta) (Griffiths, 2000). Por lo tanto, resulta más difícil determinar el sexo de un ave joven mediante métodos visuales, como durante una actividad de avistamiento de aves.

sinFalt$sex[sinFalt$sex == "Juvenil"]="Desconocido"
sinFalt$sex=droplevels(sinFalt$sex)
summary(sinFalt$sex)

Desconocido      Hembra       Macho 
        647          34          55 

Variable behavior

sinFalt$behavior=droplevels(sinFalt$behavior)
tabComp=as.data.frame(table(sinFalt$behavior))
colnames(tabComp)=c("Comportamiento","Frecuencia")
tabComp

En el caso de esta variable, se observa que hay cadenas de texto que describen más de una acción realizada por el organismo avistado, lo que ocasiona la generación de más categorías. Para este estudio, se optará por registrar en esta variable únicamente la primera acción que aparezca en la cadena de texto existente.

Se evidencia que la cadena de texto que separa la primera acción del resto del texto es ” |“, por lo que se utilizarán métodos de manipulación de cadenas de caracteres para modificar la variable.

sinFalt$behavior=as.character(sinFalt$behavior)
sinFalt$behavior=sub(" \\|.*", "", sinFalt$behavior)
sinFalt$behavior=as.factor(sinFalt$behavior)
summary(sinFalt$behavior)

               Alimentación           Autoacicalamiento 
                        103                           7 
                   Descanso Movimiento entre vegetación 
                        445                         120 
                      Vuelo 
                         61 

Con esto la variable queda reducida a 5 categorías, lo que facilitará un poco más los procesos de análisis.

Variable occurrenceRemarks

sinFalt$occurrenceRemarks=droplevels(sinFalt$occurrenceRemarks)
tabAnot=as.data.frame(table(sinFalt$occurrenceRemarks))
colnames(tabAnot)=c("Anotación","Frecuencia")
tabAnot

Esta variable debe reportar información sobre el monitoreo en el que se realizó el avistamiento, el momento del día en que se llevó a cabo, o si se avistó en un lugar particular, como un cable o un árbol muerto. Esto resulta en un alto número de categorías, con un total de 26.

Para los fines de este estudio, se registrará únicamente la información sobre el momento del día en que se realizó el avistamiento (mañana o tarde), con la intención de identificar patrones relacionados con los aspectos etológicos de las aves.

sinFalt$occurrenceRemarks=trimws(sub(".*,", "", sinFalt$occurrenceRemarks))
sinFalt$occurrenceRemarks=as.factor(sinFalt$occurrenceRemarks)
summary(sinFalt$occurrenceRemarks)

mañana  tarde 
   394    342 

Con esto, la fase de preparación se da por terminada. Se puede descargar una copia de esta versión de la base de datos en el siguiente enlace.

Análisis por categorías

Luego de la preparación de las variables, se inicia el análisis del fenómeno de biodiversidad documentado en la base de datos. Para esto, se plantean tres categorías de análisis:

1. Aspectos relacionados con la ubicación geográfica.
2. Aspectos taxonómicos de los organismos.
3. Aspectos etológicos de los organismos.

Dentro de cada una de estas categorías de análisis, se incluyen variables que se analizarán, ya sea por separado o en relación con otras, para identificar patrones y características generales del fenómeno. Una vez identificadas algunas de estas características, se procederá a realizar asociaciones entre las categorías.

Inicialmente es importante reconocer el comportamiento individual de estas variables. A continuación se analiza de manera sucinta esto, por medio de representaciones gráficas, tablas de frecuencia y medidas estadísticas (para el caso de las variables cuantitativas).

Aspectos relacionados con la ubicación geográfica.

En esta categoría se incluyen las variables decimalLongitude, decimalLatitude, parentEventID y verbatimLocality.

Variables decimalLongitude y decimalLatitude

Variable decimalLongitude

library(ggplot2)
ggplot(sinFalt, aes(x = decimalLongitude)) +
  geom_density(fill = "skyblue", color = "black", alpha = 0.7) +
  labs(title = "Longitud de ubicación de las aves",
       x = "Longitud",
       y = "Densidad") +
  theme_minimal()

summary(sinFalt$decimalLongitude)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -74.10  -74.10  -74.10  -74.09  -74.10  -74.07 

Variable decimalLatitude

ggplot(sinFalt, aes(x = decimalLatitude)) +
  geom_density(fill = "skyblue", color = "black", alpha = 0.7) +
  labs(title = "Latitud de ubicación de las aves",
       x = "Latitud",
       y = "Densidad") +
  theme_minimal()

summary(sinFalt$decimalLatitude)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.765   4.768   4.770   4.769   4.772   4.773 

Analizando estas dos variables por separado no se logra identificar comportamientos significativos en cuanto a sus distribuciones. Por ser coordenadas geográficas, tendrá sentido analizarlas en conjunto.

ggplot(sinFalt, aes(x = decimalLongitude, y = decimalLatitude)) +
  geom_point(color = "darkblue", alpha = 0.5, size = 1) +
  geom_density_2d(color = "darkgreen", alpha = 0.5) +
  labs(title = "Mapa de puntos y curvas de nivel de ubicación de aves y evento de muestreo",
       x = "Longitud",
       y = "Latitud") +
  theme_minimal()

De esto se puede ver que hay dos zonas claramente diferenciadas en las que se realizaron las actividades de observación de los organismos avistados, y esto tal vez este relacionado con la localidad. Esta relación se explorará más adelante.

Variable parentEventID

tabEvent=as.data.frame(table(sinFalt$parentEventID))
colnames(tabEvent)=c("Evento de muestreo","Frecuencia")
tabEvent

ggplot(tabEvent, aes(x = `Evento de muestreo`, y = Frecuencia)) +
  geom_bar(stat = "identity", fill = "skyblue", width = 0.5) +  
  geom_text(aes(label = Frecuencia), vjust = -0.5, size = 4, color = "black") +
  labs(title = "Evento de muestreo",
       x = "Evento",
       y = "Frecuencia") +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 8))

Variable verbatimLocality

tabLoc=as.data.frame(table(sinFalt$verbatimLocality))
colnames(tabLoc)=c("Localidad","Frecuencia")
tabLoc

ggplot(tabLoc, aes(x = "", y = Frecuencia, fill = Localidad)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  scale_fill_manual(values = c("skyblue3", "skyblue4")) +
  geom_text(aes(label = paste0(round(Frecuencia / sum(Frecuencia) * 100), "%")), 
            position = position_stack(vjust = 0.5)) +
  labs(title = "Avistamientos por Localidad") +
  theme_void() +
  theme(legend.position = "right",
        axis.text.y = element_blank(),     
        axis.ticks.y = element_blank(),    
        axis.title.y = element_blank())

Relación entre ubicación y evento de muestreo

Para la exploración de esta relación se realiza un mapa de calor bidimensional. En este mapa, se representa la densidad de aves avistadas mediante curvas de nivel, y la ubicación de los eventos de muestreo se mostrará mediante puntos, indicando dónde se identificaron las aves en dichas ubicaciones.

ggplot(sinFalt, aes(x = decimalLongitude, y = decimalLatitude)) +
  geom_density_2d(color = "red", alpha = 0.5) +
  geom_point(aes(color = parentEventID), alpha = 0.8) + 
  labs(title = "Mapa de calor de ubicación de aves y evento de muestreo",
       x = "Longitud",
       y = "Latitud") +
  theme_minimal() +
  scale_color_manual(breaks = levels(sinFalt$parentEventID), values = c("blue", "green", "purple","red"))

Se puede evidenciar que existen dos zonas de mayor concentración de los avistamientos de los organismos: el evento de muestreo E1 se concentra en una única zona (superior derecha) y los demás eventos en la otra zona de concentración (parte izquierda del gráfico).

Relación entre evento de muestreo y localidad

Para poder identificar las posibles relaciones entre estas dos variables se utiliza una tabla de contingencia, herramienta que permite relacionar dos variables cualitativas y evidenciar la distribución conjunta de las variables.

tabEventHac=table(sinFalt$parentEventID, sinFalt$verbatimLocality)
tabEventHac

     
      Hacienda la Conejera Hacienda las Mercedes
  E:1                  175                     0
  E:2                    0                   190
  E:3                    0                   189
  E:4                    0                   182

Con esto se evidencia que el evento E1 tuvo concentración en la zona de la Hacienda La Conejera, mientras que el resto de eventos se concentraron en la Hacienda Las Mercedes. También se puede ver esto en un mapa de calor, similar a la gráfica anterior.

ggplot(sinFalt, aes(x = decimalLongitude, y = decimalLatitude)) +
  geom_density_2d(color = "red", alpha = 0.5) +
  geom_point(aes(color = verbatimLocality), alpha = 0.8) + 
  labs(title = "Mapa de calor de ubicación de aves y localidad de avistamiento",
       x = "Longitud",
       y = "Latitud") +
  theme_minimal() +
  scale_color_manual(breaks = levels(sinFalt$verbatimLocality), values = c("blue", "green"))

Con esto, se puede concluir que, para facilitar el procesamiento de datos, se puede seccionar el estudio por medio de los eventos o de la localidad. Esto permitiría descartar más variables al momento de encontrar patrones y analizar relaciones entre variables de manera más sencilla. Se puede optar entonces por:

• Analizar el conjunto de datos perteneciente a la Hacienda La Conejera, lo que implica descartar las variables verbatimLocality y parentEventID, sabiendo que todos estos datos pertenecen al evento de muestreo E1.
• Analizar cada uno de los eventos restantes por separado, lo que permite descartar las variables parentEventID y verbatimLocality, teniendo en cuenta que esta última variable tomaría el valor de “Hacienda Las Mercedes” en todos los registros de estos conjunto.

Para fines de este trabajo se realizará el estudio detallado de caso de la Hacienda La Conejera.

Aspectos relacionados con la taxonomía de los organismos

En esta categoría de análisis se tendrán en cuenta las variables orden, genus, family y scientificName. Estas variables permiten reconocer la diversidad de aves observadas en el estudio. Se realizará una exploración de lo observado en toda la Reserva y lo observado en la Hacienda La Conejera. Esto puede replicarse con cada uno de los eventos de muestreo, sí es de interés del lector.

Variable order a nivel de la Reserva

sinFalt$order=droplevels(sinFalt$order)
tabOrd=as.data.frame(table(sinFalt$order))
colnames(tabOrd)=c("Orden","Frecuencia")
tabOrd = tabOrd %>%
  arrange(desc(Frecuencia))
tabOrd

tabOrd$Orden=factor(tabOrd$Orden, levels = tabOrd$Orden[order(tabOrd$Frecuencia, decreasing = FALSE)])
ggplot(tabOrd, aes(x = Frecuencia, y = `Orden`)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = Frecuencia), hjust = -0.2, size = 2, color = "black") +
  labs(title = "Ordenes de aves identificadas en la Reserva",
       x = "Frecuencia",
       y = "Ordenes") +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 8))

Variable order a nivel de la Hacienda La Conejera

AvesConejera=sinFalt[sinFalt$verbatimLocality == "Hacienda la Conejera", ]

AvesConejera$order=droplevels(AvesConejera$order)
tabOrdH=as.data.frame(table(AvesConejera$order))
colnames(tabOrdH)=c("Orden","Frecuencia")
tabOrdH = tabOrdH %>%
  arrange(desc(Frecuencia))
tabOrdH

Variable family a nivel de la Reserva

sinFalt$family=droplevels(sinFalt$family)
tabFam=as.data.frame(table(sinFalt$family))
colnames(tabFam)=c("Familia","Frecuencia")
tabFam = tabFam %>%
  arrange(desc(Frecuencia))
tabFam

tabFam$Familia=factor(tabFam$Familia, levels = tabFam$Familia[order(tabFam$Frecuencia, decreasing = FALSE)])
ggplot(tabFam, aes(x = Frecuencia, y = `Familia`)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = Frecuencia), hjust = -0.2, size = 2, color = "black") +
  labs(title = "Familias de aves identificadas en la Reserva",
       x = "Frecuencia",
       y = "Familias") +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 8))

Variable family a nivel de la Hacienda La Conejera

AvesConejera$family=droplevels(AvesConejera$family)
tabFamH=as.data.frame(table(AvesConejera$family))
colnames(tabFamH)=c("Familia","Frecuencia")
tabFamH = tabFamH %>%
  arrange(desc(Frecuencia))
tabFamH

Variable genus a nivel de la Reserva

sinFalt$genus=droplevels(sinFalt$genus)
tabGen=as.data.frame(table(sinFalt$genus))
colnames(tabGen)=c("Género","Frecuencia")
tabGen = tabGen %>%
  arrange(desc(Frecuencia))
tabGen

tabGen$Género=factor(tabGen$Género, levels = tabGen$Género[order(tabGen$Frecuencia, decreasing = FALSE)])
ggplot(tabGen, aes(x = Frecuencia, y = `Género`)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = Frecuencia), hjust = -0.2, size = 2, color = "black") +
  labs(title = "Género de aves identificadas en la Reserva",
       x = "Frecuencia",
       y = "Géneros") +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 8))

Variable genus a nivel de la Hacienda La Conejera

AvesConejera$genus=droplevels(AvesConejera$genus)
tabGenH=as.data.frame(table(AvesConejera$genus))
colnames(tabGenH)=c("Género","Frecuencia")
tabGenH = tabGenH %>%
  arrange(desc(Frecuencia))
tabGenH

Variable scientificName a nivel de la reserva

sinFalt$scientificName=droplevels(sinFalt$scientificName)
tabAv=as.data.frame(table(sinFalt$scientificName))
colnames(tabAv)=c("Especie","Frecuencia")
tabAv = tabAv %>%
  arrange(desc(Frecuencia))
tabAv

Variable scientificName a nivel de la Hacienda La Conejera

AvesConejera$scientificName=droplevels(AvesConejera$scientificName)
tabAvH=as.data.frame(table(AvesConejera$scientificName))
colnames(tabAvH)=c("Especie","Frecuencia")
tabAvH = tabAvH %>%
  arrange(desc(Frecuencia))
tabAvH

Relaciones entre variables taxonómicas

Es importante tener en cuenta que estas variables siguen una jerarquía, lo que garantiza una fuerte relación entre ellas. Por tanto, para la revisión de estas variables se utilizarán diagramas de árbol, empleando la librería treemap de R..

library(treemap)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

aves_count=sinFalt %>%
   dplyr::count(order, family, name = "Count")
treemap(aves_count,
       index = c("order", "family"),
       vSize = "Count",
       vColor = "Count",
       type = "value",
       title = "Distribución de Aves por Orden y Familia en la Reserva",
       palette = "Blues",  
       border.col = "black",  
       fontsize.labels = c(15, 12), 
       align.labels = list(c("center", "center"),c("left","top")),
       inflate.labels = FALSE)  

Por medio de esta gráfica, y para el caso de estudio, se tiene que los ordenes y familias de aves más se avistaron son:

Orden	Familia
Passeriformes	Emberizidae Turdidae Thraupidae Ictiridae Tyrannidae Fringillidae Troglodyridae Hirundinidae Parulidae
Columbiformes	Columbidae
Apodiformes	Trochilidae
Psittaciformes	Psittacidae

aves_count3=AvesConejera %>%
   dplyr::count(order, family, name = "Count")
treemap(aves_count3,
       index = c("order", "family"),
       vSize = "Count",
       vColor = "Count",
       type = "value",
       title = "Distribución de Aves por Orden y Familia en la Hacienda La Conejera",
       palette = "Blues",  
       border.col = "black",  
       fontsize.labels = c(15, 12), 
       align.labels = list(c("center", "center"),c("left","top")),
       inflate.labels = FALSE)  

Aunque en este tipo de gráfica se puede incluir una tercera variable, que en la jerarquía sería la variable genus, la representación gráfica resultante no es muy amigable para su lectura. Por ello, en este caso se opta por combinar las variables order y genus para visualizar la distribución de las aves avistadas.

aves_count2=sinFalt %>%
   dplyr::count(order, genus, name = "Count")
treemap(aves_count2,
       index = c("order", "genus"),
       vSize = "Count",
       vColor = "Count",
       type = "value",
       title = "Distribución de Aves por Orden y Género en la Reserva",
       palette = "Blues",  
       border.col = "black",  
       fontsize.labels = c(15, 12), 
       align.labels = list(c("center", "center"),c("left","top")),
       inflate.labels = FALSE)  

aves_count4=AvesConejera %>%
   dplyr::count(order, genus, name = "Count")
treemap(aves_count4,
       index = c("order", "genus"),
       vSize = "Count",
       vColor = "Count",
       type = "value",
       title = "Distribución de Aves por Orden y Género en la Hacienda La Conejera",
       palette = "Blues",  
       border.col = "black",  
       fontsize.labels = c(15, 12), 
       align.labels = list(c("center", "center"),c("left","top")),
       inflate.labels = FALSE)  

Aspectos etológicos

Esta categoría de análisis se compone de las variables occurrenceRemarks, habitat, estructuraSocial y behavior.

Variable occurrenceRemarks

tabAnot=as.data.frame(table(sinFalt$occurrenceRemarks))
colnames(tabAnot)=c("Momento del día","Frecuencia")
tabAnot

ggplot(tabAnot, aes(x = "", y = Frecuencia, fill = `Momento del día`)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  scale_fill_manual(values = c("skyblue3", "skyblue4")) +
  geom_text(aes(label = paste0(round(Frecuencia / sum(Frecuencia) * 100), "%")), 
            position = position_stack(vjust = 0.5)) +
  labs(title = "Avistamientos por momento del día-Reserva") +
  theme_void() +
  theme(legend.position = "right",
        axis.text.y = element_blank(),     
        axis.ticks.y = element_blank(),    
        axis.title.y = element_blank())

tabAnotH=as.data.frame(table(AvesConejera$occurrenceRemarks))
colnames(tabAnotH)=c("Momento del día","Frecuencia")
tabAnotH

ggplot(tabAnotH, aes(x = "", y = Frecuencia, fill = `Momento del día`)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  scale_fill_manual(values = c("skyblue3", "skyblue4")) +
  geom_text(aes(label = paste0(round(Frecuencia / sum(Frecuencia) * 100), "%")), 
            position = position_stack(vjust = 0.5)) +
  labs(title = "Avistamientos por momento del día-Hacienda La Conejera") +
  theme_void() +
  theme(legend.position = "right",
        axis.text.y = element_blank(),     
        axis.ticks.y = element_blank(),    
        axis.title.y = element_blank())

Variable habitat

tabHab=as.data.frame(table(sinFalt$habitat))
colnames(tabHab)=c("Hábitat","Frecuencia")
tabHab

ggplot(tabHab, aes(x = "", y = Frecuencia, fill =Hábitat)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  scale_fill_manual(values = c("lightblue3","skyblue3", "skyblue4","lightblue")) +
  geom_text(aes(label = paste0(round(Frecuencia / sum(Frecuencia) * 100), "%")), 
            position = position_stack(vjust = 0.5)) +
  labs(title = "Hábitat en los que se avistaron las aves-Reserva") +
  theme_void() +
  theme(legend.position = "right",
        axis.text.y = element_blank(),     
        axis.ticks.y = element_blank(),    
        axis.title.y = element_blank())

tabHabH=as.data.frame(table(AvesConejera$habitat))
colnames(tabHabH)=c("Hábitat","Frecuencia")
tabHabH

ggplot(tabHabH[tabHabH$Frecuencia > 0, ], aes(x = "", y = Frecuencia, fill =Hábitat)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0) +
  scale_fill_manual(values = c("lightblue3","skyblue3", "lightblue")) +
  geom_text(aes(label = paste0(round(Frecuencia / sum(Frecuencia) * 100), "%")), 
            position = position_stack(vjust = 0.5)) +
  labs(title = "Hábitat en los que se avistaron las aves-Hacienda La Conejera") +
  theme_void() +
  theme(legend.position = "right",
        axis.text.y = element_blank(),     
        axis.ticks.y = element_blank(),    
        axis.title.y = element_blank())

Variable estructuraSocial

tabEstr=as.data.frame(table(sinFalt$estructuraSocial))
colnames(tabEstr)=c("Est. Social","Frecuencia")
tabEstr

ggplot(tabEstr, aes(x = `Est. Social`, y = Frecuencia)) +
  geom_bar(stat = "identity", fill = "skyblue", width = 0.5) +  
  geom_text(aes(label = Frecuencia), vjust = -0.5, size = 4, color = "black") +
  labs(title = "Estructura social de las aves avistadas-Reserva",
       x = "Estr. Social",
       y = "Frecuencia") +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 8))

AvesConejera$estructuraSocial=droplevels(AvesConejera$estructuraSocial)
tabEstrH=as.data.frame(table(AvesConejera$estructuraSocial))
colnames(tabEstrH)=c("Est. Social","Frecuencia")
tabEstrH

ggplot(tabEstrH, aes(x = `Est. Social`, y = Frecuencia)) +
  geom_bar(stat = "identity", fill = "skyblue", width = 0.5) +  
  geom_text(aes(label = Frecuencia), vjust = -0.5, size = 4, color = "black") +
  labs(title = "Estructura social de las aves avistadas-Hacienda La Conejera",
       x = "Estr. Social",
       y = "Frecuencia") +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 8))

Variable behavior

sinFalt$behavior=droplevels(sinFalt$behavior)
tabComp=as.data.frame(table(sinFalt$behavior))
colnames(tabComp)=c("Comportamiento","Frecuencia")
tabComp

ggplot(tabComp, aes(x = Comportamiento, y = Frecuencia)) +
  geom_bar(stat = "identity", fill = "skyblue", width = 0.5) +  
  geom_text(aes(label = Frecuencia), vjust = -0.5, size = 4, color = "black") +
  labs(title = "Comportamiento de las aves avistadas-Reserva",
       x = "Comportamiento",
       y = "Frecuencia") +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 8))

AvesConejera$behavior=droplevels(AvesConejera$behavior)
tabCompH=as.data.frame(table(AvesConejera$behavior))
colnames(tabCompH)=c("Comportamiento","Frecuencia")
tabCompH

ggplot(tabCompH, aes(x = Comportamiento, y = Frecuencia)) +
  geom_bar(stat = "identity", fill = "skyblue", width = 0.5) +  
  geom_text(aes(label = Frecuencia), vjust = -0.5, size = 4, color = "black") +
  labs(title = "Comportamiento de las aves avistadas-Hacienda La Conejera",
       x = "Comportamiento",
       y = "Frecuencia") +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 8))

Relaciones entre las variables etológicas

Para examinar cómo interactúan estas variables, se realizará un diagrama de árbol para explorar la partición correspondiente a la Hacienda La Conejera en los avistamientos realizados en la mañana, aprovechando su potencial para resumir el comportamiento de las variables en conjunto. A diferencia de las variables taxonómicas, estas no siguen una jerarquía específica, por lo que queda a discreción del investigador elegir el orden en que se presentarán en el diagrama. Para el presente estudio, se trabajarán estas variables en el siguiente orden: comportamiento, hábitat y estructura social.

AvesConeAM=AvesConejera[AvesConejera$occurrenceRemarks == "mañana", ]

AvesConeAM_count=AvesConeAM %>%
  dplyr::count(behavior,habitat,estructuraSocial,name = "Count")
treemap(AvesConeAM_count,
        index = c( "behavior","habitat","estructuraSocial"),
        vSize = "Count",
         vColor = "Count",
         type = "value",
         title = "Distribución de las aves según su comportamiento, hábitat y estr. social",
         palette = "Blues",  
         border.col = "black",  
         fontsize.labels = c(16, 13, 10), 
         align.labels =list(
           c("center", "top"),   
           c("left", "center"),  
           c("right", "bottom")  
         ),
         inflate.labels = FALSE)

Como se puede ver, la mayoría de los avistamientos de la mañana en la Hacienda La Conejera correspondieron a aves descansando en hábitats de bosque mixto nativo, en grupos coespecíficos. Posteriormente, se observó una cantidad considerable de aves desplazándose entre la vegetación, principalmente en hábitats de bosque nativo, en grupos coespecíficos o en parejas. En menor proporción, se avistaron aves alimentándose, principalmente en bosque nativo y en parejas. Con menor frecuencia, se observaron aves volando, en igual proporción en bosques mixtos y nativos, usualmente en grupos coespecíficos.

Análisis de asociaciónes por categorías

Una vez estudiado el fenómeno de biodiversidad por categorías, se propone buscar posibles relaciones entre estas para extraer información adicional que permita identificar patrones en el muestreo recolectado en la base de datos. Para esto, se utilizará un método estadístico multivariado conocido como Análisis de Correspondencia Múltiple (ACM). En el presente estudio, se analizará la partición correspondiente a la localidad de la Hacienda La Conejera y las observaciones realizadas en la jornada de la mañana. Las variables consideradas serán el orden (order), el comportamiento (behavior), la estructura social (estructuraSocial) y el hábitat (habitat).

Análisis de correspondencia múltiple

library(FactoMineR)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

datosACM=subset(AvesConeAM, select = c("order",
                                    "behavior",
                                    "estructuraSocial",
                                    "habitat"))
mca=MCA(datosACM,graph=FALSE)

library(factoextra)

Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

library(pander)
eigenval=get_eigenvalue(mca)
pander(eigenval)

--------------------------------------------------------------------------
         eigenvalue   variance.percent   cumulative.variance.percent 
------------ ------------ ------------------ -----------------------------
 **Dim.1**      0.5181          20.72                    20.72            

 **Dim.2**      0.3563          14.25                    34.97            

 **Dim.3**      0.3418          13.67                    48.65            

 **Dim.4**      0.2924           11.7                    60.34            

 **Dim.5**      0.2597          10.39                    70.73            

 **Dim.6**      0.1932          7.728                    78.46            

 **Dim.7**      0.1764          7.056                    85.51            

 **Dim.8**      0.1599          6.397                    91.91            

 **Dim.9**      0.1075          4.301                    96.21            

 **Dim.10**    0.09471          3.789                     100             
--------------------------------------------------------------------------

scpl=fviz_screeplot(mca,addlabels=TRUE)
scpl

Se tiene que las primeras cinco dimensiones determinadas por el ACM explican un 70.73% de la varianza del fénomeno estudiado. Se procede a explorar las actegorías que definen estas dimensiones.

varCont=mca$var$contrib
varContT=rowSums(varCont[, 1:5])
varContP=100*varContT/sum(varContT)
varContDF=data.frame(Variable=names(varContP), Contribución=varContP)
varContDF=varContDF[order(varContDF$Contribución, decreasing = TRUE), ]
varContDF

par(mar = c(15, 4, 4, 2) + 0.1)
barplot=barplot(varContDF$Contribución, 
                         names.arg = varContDF$Variable, 
                         las = 2, 
                         col = "#4682b4",
                         main = "Contribución de las categorias de las variables a las Dimensiones 1-5 (%)",
                         ylim = c(0, max(varContDF$Contribución) * 1.2))
text(x = barplot,
      y = varContDF$Contribución,
      label = sprintf("%.1f%%", varContDF$Contribución),
      pos = 3, cex = 0.8, col = "black")

Al detallar las categorías que determinan estas cinco dimensiones, se encuentran que entre las principales 14 se resaltan ordenes taxonómicos como los cuculiformes, galliformes y los apodiformes, así como comportamientos como volar, alimentarse y moverse entre la vegetación, y la estancia de las aves en parejas. Además, hábitats como el bosque mixto nativo y el bosque nativo también figuran entre las categorías más significativas

Clusterización por k-means

Aprovechando el proceso realizado con el ACM, se procede a hacer un proceso de clusterización por k-means, con el fin de identificar grupos con características similares. Para este proceso, se hará uso de las primeras dos dimensiones generadas por el ACM (que explican el 35% de la varianza). Inicialmente se revisa las categorías que determinan estas dos dimensiones.

varContT2=rowSums(varCont[, 1:2])
varContP2=100*varContT2/sum(varContT2)
varContDF2=data.frame(Variable=names(varContP2), Contribución=varContP2)
varContDF2=varContDF2[order(varContDF2$Contribución, decreasing = TRUE), ]
varContDF2

par(mar = c(15, 4, 4, 2) + 0.1)
barplot2 = barplot(varContDF2$Contribución, 
                   names.arg = varContDF2$Variable, 
                   las = 2, 
                   col = "#4682b4",
                   main = "Contribución de las categorías de las variables a las Dimensiones 1-2 (%)",
                   ylim = c(0, max(varContDF2$Contribución) * 1.2))
text(x = barplot2,
     y = varContDF2$Contribución,
     label = sprintf("%.1f%%", varContDF2$Contribución),
     pos = 3, cex = 0.8, col = "black")

Entre las categorías que tienen contribuciones significativas a estas dimensiones se encuentran los órdenes Galliformes, Columbiformes y Cuculiformes, así como comportamientos como la alimentación y el movimiento entre vegetación, la estancia en grupos coespecíficos o en solitario.

Luego de esto, se procede a la exploración de asociaciones entre observaciones y categorías por medio de un diagrama de biplot.

par(mar = c(15, 4, 4, 2) + 0.1)
barplot2 = barplot(varContDF2$Contribución, 
                   names.arg = varContDF2$Variable, 
                   las = 2, 
                   col = "#4682b4",
                   main = "Contribución de las categorías de las variables a las Dimensiones 1-2 (%)",
                   ylim = c(0, max(varContDF2$Contribución) * 1.2))
text(x = barplot2,
     y = varContDF2$Contribución,
     label = sprintf("%.1f%%", varContDF2$Contribución),
     pos = 3, cex = 0.8, col = "black")

![](dm12nb)

Error: inesperado '[' en "!["

Usando este gráfico como base, se inicia el proceso de clusterización.

library(cluster)
indCoords=mca$ind$coord
set.seed(123)
kmRes=kmeans(indCoords, centers = 3, nstart = 25)
mca$ind$cluster=kmRes$cluster
indData=data.frame(indCoords,cluster = as.factor(kmRes$cluster))

createPolygon=function(cluster_data){
  hull=chull(cluster_data$Dim.1, cluster_data$Dim.2)
  hull=c(hull, hull[1])
  cluster_data[hull, ]
}

polygonData=indData %>%
  group_by(cluster) %>%
  do(createPolygon(.))

ggplot(indData, aes(x = Dim.1, y = Dim.2, color = cluster)) +
  geom_point() +
  geom_polygon(data = polygonData, aes(x = Dim.1, y = Dim.2, group = cluster, fill = cluster), alpha = 0.2) +
  theme_minimal() +
  labs(title = "Clusterización de Aves", x = "Dimensión 1", y = "Dimensión 2")

combPlot +
  geom_point(data = indData, aes(x = Dim.1, y = Dim.2, color = cluster)) +
  geom_polygon(data = polygonData, aes(x = Dim.1, y = Dim.2, group = cluster, fill = cluster), alpha = 0.2)

Este proceso se realizó estableciendo tres clústeres, de los cuales dos están claramente definidos en la gráfica. El primer grupo muestra que las aves del orden Columbiformes tienden a habitar bosques mixtos nativos y son avistadas en grupos coespecíficos. En el segundo clúster, se observa que las aves de los órdenes Passeriformes y Apodiformes suelen habitar bosques nativos y son avistadas mientras descansan en parejas o solas. El último grupo no está tan bien delimitado como los dos anteriores, mostrando una alta heterogeneidad y no agrupándose según las categorías predominantes de las dos dimensiones analizadas en este proceso.

Referencias

Barrera, M. C. A. (2023, agosto). Caracterización de avifauna - Reserva Thomas Vander Hammen v2.1. Dataset/Occurrence. https://doi.org/10.15472/alrlwe

Griffiths, Richard, and D Phil. (2000). “Sex Identification in Birds.”