Introducción

La ecología busca comprender las interacciones entre los organismos y el ambiente, así como identificar las reglas generales que dilucidan el funcionamiento de los ecosistemas, para en última instancia mejorar nuestra capacidad para predecir los cambios que se producen en ellos. En Ecología Marina, tanto para los hábitats pelágicos como para los bentónicos, las preguntas importantes siguen siendo las mismas: (1) ¿Cuáles son los procesos que controlan la estructura y el funcionamiento de los ecosistemas marinos? (2) ¿Qué patrones ecológicos surgen a distintas escalas espacio-temporales y cuáles son sus principales impulsores? (3) ¿Cómo responderán los organismos marinos a las crecientes presiones antropogénicas?. Recientemente se resaltó el potencial de la ecología basada en rasgos para el estudio de los ecosistemas acuáticos. Los enfoques basados en los rasgos (definidos en la investigación ecológica como cualquier método que se centra en las propiedades individuales de los organismos) en lugar de las especies podrían proporcionar un marco común para abordar las preguntas formuladas arriba. Esta guía de ejemplos y ejercicios en R proporciona los conocimientos técnicos iniciales a los análisis de datos reales de rasgos en combinación con datos comunitarios, ambientales o filogenéticos.

Datos de rasgos

Paquetes necesarios y funciones adicionales


# install.packages("devtools")
library(devtools)
install_github("larsito/traitor")

# install.packages("taxize")
library(taxize)
# install.packages("Taxonstand")
library(Taxonstand)

# install.packages("traits")
library(traits)

# install.packages("plyr")
library(plyr)

Obtención de datos de rasgos (datos LEDA)

Hemos visto antes que existen muchos repositorios en internet de los que podemos obtener datos. Algunas bases de datos pueden estar ubicadas físicamente en un servidor al que no tenemos acceso directo, y que requieren hacer una petición a una institución que gestione los datos. En el ejemplo que utilizamos a continuación, los datos son accesibles de forma más directa, en forma de archivo .txt descargable desde un sitio web dedicado a ello. Sea cual sea la fuente, la utilización de datos procedentes de bases de datos o conjuntos de datos de artículos publicados suele requerir cierta “limpieza” necesaria, como veremos.

Para comenzar utilizaremos un único rasgo de la base de datos LEDA. Descargaremos, a través de R, un único archivo .txt que contiene datos sobre la altura del dosel, es decir, la altura sobre la superficie del suelo de las hojas fotosintéticamente activas más altas. Utilizando la función “download.file”, hay dos rutas que debemos especificar. La primera, a través del argumento “url”, es la ubicación del archivo de texto en Internet (es decir, una dirección de enlace). El argumento “destfile” se utiliza para especificar la ruta a la que descargar el archivo, incluyendo el nombre del mismo. Si sólo especificamos el nombre del archivo .txt, el destino será automáticamente el directorio de trabajo actual (que es el predeterminado cuando no se especifica ninguna otra ruta antes del nombre del archivo).

La función proporciona información e incluso una barra de progreso que indica la parte del archivo que se ha descargado. En caso de una conexión a Internet normal, esta descarga debería ser bastante rápida y realizarse en unos pocos segundos.


urlLeda <- 
  "https://uol.de/f/5/inst/biologie/ag/landeco/download/LEDA/Data_files/canopy_height.txt"
miurl <- "./0Data/EcolMarina/height_leda.txt"

download.file(url = urlLeda, destfile = miurl)

Lo primero que hay que hacer es abrir el archivo en un editor de texto de su elección y echar un vistazo al contenido. A continuación, utilicen read.table con el argumento “skip” fijado en 3 (el número de líneas que deben saltarse, aunque en su editor puedan parecer más de 3 líneas). La tabla importada debe guardarse como un objeto, aquí llamado “height_leda”. También podemos ver en el editor que se utilizan puntos y comas para separar las columnas y tenemos una fila que contiene nombres de columnas, por lo que establecemos los argumentos “sep” y “header” en consecuencia. Además, nos deshacemos de todos los llamados espacios en blanco.


height_leda <- read.table(miurl,
                          skip = 3, sep = ";",
                          header = TRUE,
                          strip.white = TRUE)

# View(height_leda) # visualiza el data frame completo en una pestaña aparte

La función “read.table” cambia por defecto todas las celdas vacías en valores NA. La función read.table tiene aproximadamente dos docenas de argumentos que pueden utilizarse para ajustar la forma en que se importan los datos. Recomiendo que la estudien por si alguna vez se encuentran con algún problema específico al importar los datos.
La base especifica, entre otras cosas, información sobre el valor medio del rasgo, los valores mínimos y máximos, el tamaño de la muestra en que se basan estas medidas, las referencias (en caso de que los datos se hayan extraído de la literatura publicada o de bases de datos ya disponibles), y la metodología en que se basa la medición. Potencialmente, toda esta información adicional puede utilizarse para optimizar nuestros datos. Por ejemplo, cuando los datos obtenidos de “experimentos de laboratorio/invernadero/jardín” (véase la columna ‘método.general’) no deben incluirse porque estas plantas podrían no haber sido cultivadas en condiciones “naturales”.

La base de datos de LEDA deja claro lo importante que es tener algún tipo de clave que nos ayude a entender los datos. Como mínimo, deberíamos tener algunos metadatos (datos sobre los datos) que nos informen sobre las abreviaturas utilizadas, las metodologías aplicadas, etc. A menudo, estos metadatos se almacenan como archivos adicionales .txt o Excel que acompañan a los archivos de datos de rasgos reales. LEDA cuenta con una descripción muy detallada de cómo se midieron los diferentes rasgos y cómo se representan en la base de datos. Noten que para muchas especies hay más de un registro de altura en LEDA. Esto suele deberse a que la información para la misma especie se encontró en diferentes fuentes, y todas ellas se introdujeron en la base de datos. Si estamos interesados en obtener un único valor por especie, necesitamos agregar los datos a nivel de especie. Como queremos el valor medio de la altura de una especie, utilizamos la media.


height_leda_agg <- aggregate(height_leda$single.value..m., # la variable u
                             # objeto que desea agregar
                             list(height_leda$SBS.name), # nivel de
                             # agregación (usualmente categórica)
                             FUN = mean)
head(height_leda_agg)

Bases con información taxonómica

Para algunos grupos taxonómicos, existe la posibilidad de acceder a bases de datos taxonómicas que sirven de referencia estándar para corroborar que las especies de la lista estén nombradas adecuadamente y se correspondan con los rasgos descritos. Por lo general, estas bases de datos ofrecen la posibilidad de buscar el nombre de una especie concreta y recuperar información relevante sobre su situación y sus sinónimos. Varios paquetes de R ofrecen la posibilidad de consultar estos sitios web directamente desde R.
El paquete taxize permite acceder a un gran número de bases de datos de taxones, y ofrece una gran flexibilidad en el tipo de información que se puede recuperar de esas bases. El número de funciones es un poco abrumador al principio, pero son básicamente una combinación de dos partes. La primera parte es una abreviatura de la base de datos de taxones que se quiere consultar, y la segunda, el conjunto de funciones específicas que se utilizan para obtener información de cada base de datos. El tipo de información que se puede extraer difiere entre las bases de datos. Es posible que queramos ver si realmente estamos utilizando nombres de especies que estén de acuerdo con una determinada fuente de datos taxonómicos. Esto se puede hacer con la función “resolve”.


resolve("Vulpes vulpes")
#> $gnr
#> [1] "Error: no data found"

Así, pueden estar bastante seguras de que hemos buscado un nombre de especie que se utiliza realmente como tal en una amplia variedad de fuentes. No obstante, no obtenemos mucha más información que esa (aparte de saber en qué bases de datos se utiliza un nombre determinado). Algunas fuentes de datos taxonómicos que están disponibles a través de taxize contienen información sobre grupos taxonómicos específicos. Por ejemplo, al trabajar con organismos marinos podemos utilizar taxize para consultar el “Registro Mundial de Especies Marinas” (WoRMS).

Nota al margen: Comenté los “chunks” correspondientes al uso de Worms porque LaTex en RMarkdown no puede manejar el símbolo de “visto bueno” de la salida del primer chunk.


# mono.mono_ID <- get_wormsid("Monodon monoceros")
# mono.mono_ID

Este número de identificación se utiliza para llamar a funciones adicionales que trabajan en la base de datos WoRMS. Por ejemplo, para obtener una visión general de la clasificación taxonómica de las especies a partir del reino.


# classification(mono.mono_ID)

O para obtener un data.frame que muestre todos los diferentes sinónimos de esa especie, incluyendo la autoría de estos sinónimos, información sobre la clasificación, etc.


# synonyms(mono.mono_ID)

La información que dará la función taxize depende de la base de datos que se consulte y no ofrece las soluciones más sencillas para resolver los problemas de sinonimia. Un paquete útil a este respecto es Taxonstand, aunque tiene la limitación de que sólo se ocupa de las especies vegetales. Taxonstand también ofrece una interfaz con una base de datos de taxones, concretamente The Plant List. La función “TPL” puede buscar una lista completa de especies, mientras que “TPLck” sólo busca da una especie. TPL puede manejar hasta cierto punto pequeños errores tipográficos en el nombre específico, pero no en el nombre del género. Una característica útil de Taxonstand es que proporciona una salida detallada de las especies sobre su estado taxonómico en forma tabulada:


Taxonstand::TPL(c(sp = "Elymus repens"))

Cotejando las especies en matrices de rasgos y comunitarias

El uso de “Taxonstand” y su función “TPL” para traducir nuestra lista de especies a nombres de especies aceptados es muy útil. Nos permite cotejar diferentes data.tables basados en especies, incluso cuando las diferentes data.tables contienen la misma especie pero con un nombre diferente o con (un grado razonable de) errores. Lo ejemplificamos construyendo un conjunto de datos artificial. En primer lugar, creamos una tabla de datos de rasgos (“height”) con cuatro especies de plantas, con nombres que no son todos correctos. Los nombres correctos de las especies utilizadas en el ejemplo son “Abies alba” (error tipográfico), “Bellis dubia” (sinónimo), “Helianthus annua” (error ortográfico), mientras que “Cocos nucifera” no está mal escrito ni es un sinónimo. A continuación, construimos los datos de la comunidad de presencia/ausencia, que deberían incluir las mismas especies, pero que, comparados con nuestra primera tabla de rasgos, tienen todos los nombres correctos (“com”).


height <- data.frame(species = c("Abies alpa", "Chrysanthemum bellidiflorum",
                                 "Helianthus annua", "Cocos nucifera"),
                     height = c(45, 0.1, 1.9, 28))
height


com <- matrix(c(1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0), 4, 4)
colnames(com) <- c("Abies alba", "Bellis dubia", "Helianthus annuus",
                   "Cocos nucifera")
rownames(com) <- c("community 1", "community 2", "community 3",
                   "community 4")
com
#>             Abies alba Bellis dubia Helianthus annuus Cocos nucifera
#> community 1          1            0                 0              1
#> community 2          1            1                 0              0
#> community 3          0            1                 1              1
#> community 4          0            1                 1              0

El primer paso para cotejar los datos de la comunidad con los datos de los rasgos es ver si los nombres de las especies en ambos conjuntos de datos son los mismos. Se pueden utilizar dos funciones útiles en R, “intersect” y “setdiff”. “intersect” nos da la intersección de dos variables (por ejemplo, la intersección de las variables $\{a, b, c\}$ y $\{c, d, e\}$ sería $\{c\}$). Como queremos obtener el conjunto de nombres de especies comunes en los dos conjuntos de datos (comunidad y rasgos), podemos utilizar la intersección en los nombres de las columnas de los datos de la comunidad y de las especies de los datos de los rasgos. La función “setdiff” resulta útil para detectar qué especies están presentes en un conjunto de datos, pero no en el otro (ojo, “setdiff” devuelve las especies que están presentes en el primer objeto, pero no en el segundo - si se alterna el orden, el resultado es diferente). Si encontramos inconcistencias en los nombres de las especies entre ambos conjuntos de datos, enviamos una consulta con las listas de nombres de especies sospechosas usando la función “TPL”.


intersect(colnames(com), height$species)
#> [1] "Cocos nucifera"
setdiff(colnames(com), height$species)
#> [1] "Abies alba"        "Bellis dubia"      "Helianthus annuus"

tpl_com <- TPL(colnames(com))
tpl_com


tpl_height <- TPL(height$species)
tpl_height

Para poder vincular los dos conjuntos de datos, los nombres defectuosos de las especies deben sustituirse por los nombres adecuados. Para ello, puede utilizarse una combinación de las variables New.Genus y New.Species. Noten que para obtener los nombres de las especies tenemos que combinar (es decir, usar paste) los nombres de los géneros y las especies. Luego, pueden añadirse, o sustituir por, los nombres adecuados. Hecho eso utilizamos una vez más la función de intersección y verificamos que todos los nombres de las especies coincidan.


new_species_names_com <- paste(tpl_com$New.Genus, tpl_com$New.Species)
new_species_names_height <- paste(tpl_height$New.Genus, tpl_height$New.Species)

height$new_species_names_height <- new_species_names_height
colnames(com) <- new_species_names_com
intersect(colnames(com), height$new_species_names_height)
#> [1] "Abies alba"        "Bellis dubia"      "Helianthus annuus"
#> [4] "Cocos nucifera"

Uso del paquete “traitor” para apoyar su campaña de muestreo de rasgos

Los datos faltantes pueden afectar al resultado de los análisis de rasgos, especialmente cuando los datos de rasgos faltan para especies que son abundantes en una comunidad. El número de “NA” que podemos permitirnos en nuestros datos, y seguir considerando que los resultados del análisis son correctos, depende del contexto. Sin embargo, una regla general sugerida es que no sea más del 20% de la abundancia o biomasa total. Una estrategia para mitigar el efecto de los datos que faltan es cuantificar la cantidad de información que se pierde, y después medir los rasgos de determinadas especies. El paquete “traitor” puede utilizarse para priorizar qué especies deben medirse. Comenzamos cargando el paquete y un conjunto de datos asociado.


library(traitor)
data(ohrazeni)

La función “sampleSpecies” puede usarse para identificar de cuáles especies debemos obtener información de sus rasgos, asumiendo que tenemos datos faltantes. La función tiene cuatro argumentos: el argumento “com” proporciona los datos de la comunidad; mientras que “avail_sp” nos dice para qué especies de la comunidad contamos con sufieciente información de sus rasgos (y para cuáles no). Basándose en esta información “sampleSpecies” sugiere qué especies de la comunidad debemos muestrear. El argumento “thresh” define el umbral mínimo de abundancia sumada de las especies para las que necesitamos datos de rasgos. Por ejemplo, fijando este umbral en el 80% de la abundancia total, y suponiendo que en una parcela dada hay datos de rasgos para las especies que constituyen el 70% de la abundancia relativa, necesitamos muestrear los rasgos de tantas especies como sea necesario para un 10% adicional de la abundancia total. Puede tratarse de una sola especie o de más de una, en función de su abundancia. Para entender el cuarto argumento, secuencial, tenemos que introducir los escenarios “pool-wise” (Escenario 1) y “plot-wise” (Escenario 2), que son aplicados automáticamente por “sampleSpecies”.
Imaginemos una matriz de comunidad de tres parcelas. Con la información de “avail_sp” podemos calcular para cada parcela qué especies deben ser muestreadas para alcanzar el umbral de abundancia total del 80%, definido en el argumento “thresh”. Dado que la composición de especies será diferente entre las parcelas, “sampleSpecies” podría sugerir diferentes especies para alcanzar el 80% en cada una de las parcelas. Veamos esto con unos datos simples a modo de ejemplo. Por el momento, fijamos el argumento “sequential” en “FALSE” (volveremos sobre su significado más adelante).


# 1. Crear una comunidad ficticia con 3 muestras (filas) y 5 especies (columnas)
com_toy <- matrix(c(50, 40, 30, 20, 0, 30, 0, 10, 0, 5, 10, 0, 0, 0, 5), 3, 5)

# 2. Dar los "nombres" sp1 a sp5 a las especies
colnames(com_toy) <- as.character(paste("sp", c(1:5)))

# 3. Crear un vector que defina para qué especies hay datos
# de rasgos (1) y para las que no los hay (0)
avail <- c(1, 0, 1, 0, 1)

# 4. Ejecutar la función sampleSpecies en los datos ficticios
sampleSpecies(com = com_toy, avail_sp = avail, thresh = 0.8, sequential = FALSE)
#> Warning: avail_sp is without names. Species are assumed to be in the same order
#> as in community data.
#> $Original
#> $Original$com_available
#>      sp 1 sp 3 sp 5
#> [1,]   50    0    0
#> [2,]   40   10    0
#> [3,]   30    0    5
#> 
#> $Original$available_pool
#> [1] 0.675
#> 
#> $Original$available_plot
#> [1] 0.6666667 0.8333333 0.5384615
#> 
#> 
#> $Scenario1
#> $Scenario1$sample_species
#> [1] "sp 2"
#> 
#> $Scenario1$com_available
#>      sp 1 sp 2 sp 3 sp 4 sp 5
#> [1,]   50   20    0   NA    0
#> [2,]   40    0   10   NA    0
#> [3,]   30   30    0   NA    5
#> 
#> $Scenario1$available_pool
#> [1] 0.925
#> 
#> $Scenario1$available_plot
#> [1] 0.9333333 0.8333333 1.0000000
#> 
#> 
#> $Scenario2
#> $Scenario2$sample_species
#> $Scenario2$sample_species[[1]]
#> [1] "sp 2"
#> 
#> $Scenario2$sample_species[[2]]
#> [1] "none"
#> 
#> $Scenario2$sample_species[[3]]
#> [1] "sp 2"
#> 
#> 
#> $Scenario2$com_available
#>      sp 1 sp 2 sp 3 sp 4 sp 5
#> [1,]   50   20    0   NA    0
#> [2,]   40   NA   10   NA    0
#> [3,]   30   30    0   NA    5
#> 
#> $Scenario2$available_pool
#> [1] 0.925
#> 
#> $Scenario2$available_plot
#> [1] 0.9333333 0.8333333 1.0000000

El resultado de la función se almacena en una lista con tres elementos, llamados Original, Escenario1 y Escenario2, respectivamente. El elemento Original contiene una matriz de datos de la comunidad, que muestra sólo las especies para las que hay datos de rasgos, de modo que en nuestro caso sólo tiene tres especies (sp1, sp3 y sp5) en lugar de las cinco especies. Las especies 1, 3 y 5 representan conjuntamente el 67,5 % de la abundancia total en el estudio global (véase “Original.available_pool”). Sólo la segunda parcela cumpliría el umbral de abundancia total del 80% sin el muestreo de especies adicionales (véase “Original.available_plot”).

Los elementos Escenario1 y Escenario2, cada uno con una entrada adicional, especifican qué especies debemos muestrear en el escenario dado. Para alcanzar el umbral de abundancia del 80% en el escenario por grupos, se debe muestrear la especie 2. Si es así, con esta especie en los datos se cubre el 92,5% de la abundancia total (véase “Scenario1.available_pool”). Esto es más de lo necesario porque la función selecciona la especie más abundante que ayuda a alcanzar el valor del umbral, maximizando el beneficio que obtenemos del muestreo de la menor cantidad de especies adicionales.

En la matriz de datos de la comunidad bajo “.com_available” vemos ahora la comunidad representada como si ya hubiéramos muestreado esa especie número 2, es decir, podemos utilizar sus datos de abundancia en las parcelas para calcular la diversidad funcional, etc. Obsérvese también que, aunque todavía no hemos considerado las parcelas individuales para averiguar qué especies muestrear (es decir, el Escenario 2), ya podemos ver que también las abundancias de las especies a nivel de parcela (“.Escenario1.disponible_parcela”) han superado todas el 80%, por lo que hemos matado dos pájaros de un tiro al utilizar sólo el Escenario 1. Esta es también la razón por la que la salida en el tercer elemento “.Scenario2” es idéntica a la del Escenario 1, con algunas pequeñas excepciones. Para el escenario por parcelas, obtenemos sugerencias para cada parcela individual sobre qué especies muestrear. La función sugiere medir la especie 2 para las parcelas 1 y 3, pero “ninguna” para la muestra 2. En la muestra 2, el 83% de la abundancia total ya está representada por especies con datos de rasgos.
Procedemos a realizar un pequeño cambio en el objeto “avail” que hemos creado, para ver las consecuencias al plantear diferencias importantes entre el Escenario 1 y el Escenario 2, mediante la inclusión de una especie adicional son datos de rasgos.


# Crear un vector que defina los datos de los rasgos que faltan para una
# especie adicional
avail2 <- c(1, 0, 0, 0, 1)
# Volver a ejecutar la función sampleSpecies en los datos ficticios
sampleSpecies(com = com_toy, avail_sp = avail2, thresh = 0.8,
              sequential = FALSE)
#> Warning: avail_sp is without names. Species are assumed to be in the same order
#> as in community data.
#> $Original
#> $Original$com_available
#>      sp 1 sp 5
#> [1,]   50    0
#> [2,]   40    0
#> [3,]   30    5
#> 
#> $Original$available_pool
#> [1] 0.625
#> 
#> $Original$available_plot
#> [1] 0.6666667 0.6666667 0.5384615
#> 
#> 
#> $Scenario1
#> $Scenario1$sample_species
#> [1] "sp 2"
#> 
#> $Scenario1$com_available
#>      sp 1 sp 2 sp 3 sp 4 sp 5
#> [1,]   50   20   NA   NA    0
#> [2,]   40    0   NA   NA    0
#> [3,]   30   30   NA   NA    5
#> 
#> $Scenario1$available_pool
#> [1] 0.875
#> 
#> $Scenario1$available_plot
#> [1] 0.9333333 0.6666667 1.0000000
#> 
#> 
#> $Scenario2
#> $Scenario2$sample_species
#> $Scenario2$sample_species[[1]]
#> [1] "sp 2"
#> 
#> $Scenario2$sample_species[[2]]
#> [1] "sp 4"
#> 
#> $Scenario2$sample_species[[3]]
#> [1] "sp 2"
#> 
#> 
#> $Scenario2$com_available
#>      sp 1 sp 2 sp 3 sp 4 sp 5
#> [1,]   50   20   NA   NA    0
#> [2,]   40   NA   NA   10    0
#> [3,]   30   30   NA   NA    5
#> 
#> $Scenario2$available_pool
#> [1] 0.925
#> 
#> $Scenario2$available_plot
#> [1] 0.9333333 0.8333333 1.0000000

Ahora nos encontramos con un 63% de la abundancia relativa cubierta por especies con datos de rasgos a nivel de pool, y casi con sólo la mitad (54%) en la tercera parcela. Ciertamente, no sería lo más fiable calcular ningún índice de diversidad o composición basado en rasgos a partir de datos que representan una proporción tan pequeña de la vegetación muestreada. Con el escenario de pool-wise , podemos remediar esta situación, como antes, muestreando rasgos para la especie 2, que cumple nuestro umbral del 80%, pero con una abundancia relativa total menor que en el primer ejemplo (87,5% frente a 92,5%). En este escenario, la segunda parcela seguiría teniendo sólo datos de las especies que constituyen dos tercios de la abundancia relativa de la parcela. Aquí es donde entra en juego el Escenario 2. Para la segunda parcela, sugiere muestrear la especie 4 en lugar de la especie 2, que es la especie que se debe muestrear para las parcelas 1 y 3. Esto se debe a que la especie 2 no está presente en absoluto en la muestra 2, por lo que no nos ayuda a alcanzar el umbral para esa parcela. Por otro lado, la especie 4 representa el 17% de la abundancia total de 60 en esa parcela, y sumando eso al 67% o abundancia relativa de la especie 1 alcanzamos el umbral del 80%.

Por último, pasemos al argumento “sequential”. Cuando el argumento “sequential” se establece como TRUE, la función determina primero qué especies deben ser muestreadas para alcanzar el umbral según el Escenario 1, es decir, qué especies deben ser muestreadas para garantizar que la abundancia global (agrupada) cubierta sea al menos el valor establecido por “thresh”. Sin embargo, dado que esto lleva a que las parcelas individuales no necesariamente alcancen este umbral (como hemos visto anteriormente), posteriormente se determinan las especies por el Escenario 2 (es decir, por parcela) que necesitan ser muestreadas además de las determinadas por el Escenario 1, con el fin de alcanzar el umbral para cada una de las parcelas. Esta estrategia de muestreo sería adecuada si se necesitan valores de rasgos comunitarios (no importando de qué parcela individual procedan las mediciones de rasgos). De hecho, podrían ser incluso de individuos medidos fuera de la muestra. También puede ocurrir que el umbral ya se haya alcanzado para cada parcela con el Escenario 1 (como en nuestro primer ejemplo), en cuyo caso se muestra una advertencia y el valor para el Escenario 2 en la salida es NULL, pero sólo si el argumento “sequential” está efectivamente establecido en TRUE.


sampleSpecies(com = com_toy, avail_sp = avail, thresh = 0.8,
              sequential = TRUE)
#> Warning: avail_sp is without names. Species are assumed to be in the same order
#> as in community data.
#> Warning: Threshold reached by Scenario 1 alone
#> $Original
#> $Original$com_available
#>      sp 1 sp 3 sp 5
#> [1,]   50    0    0
#> [2,]   40   10    0
#> [3,]   30    0    5
#> 
#> $Original$available_pool
#> [1] 0.675
#> 
#> $Original$available_plot
#> [1] 0.6666667 0.8333333 0.5384615
#> 
#> 
#> $Scenario1
#> $Scenario1$sample_species
#> [1] "sp 2"
#> 
#> $Scenario1$com_available
#>      sp 1 sp 2 sp 3 sp 4 sp 5
#> [1,]   50   20    0   NA    0
#> [2,]   40    0   10   NA    0
#> [3,]   30   30    0   NA    5
#> 
#> $Scenario1$available_pool
#> [1] 0.925
#> 
#> $Scenario1$available_plot
#> [1] 0.9333333 0.8333333 1.0000000

Si el sequential se fijara en FALSE, las especies a muestrear según el Escenario 2 son independientes del Escenario 1. Las especies indicadas para ser medidas por el Escenario 1 son las mismas que cuando el secuencial se establece en TRUE, pero el Escenario 2 indica las especies sin tener en cuenta el Escenario 1. Por lo tanto, la función proporciona ahora un conjunto completo de especies que deben ser muestreadas en cada parcela individual del estudio. Esta sería una estrategia de muestreo adecuada cuando, por ejemplo, se quisiera tener mediciones de rasgos por parcela, es decir, se podría tener en cuenta la variación de rasgos en el nivel de la parcela. Es importante señalar que estas funciones necesitan algún tipo de datos de abundancia (semi)cuantitativos (por ejemplo, frecuencia, valores de cobertura, recuentos individuales, biomasa) y no funcionan con matrices comunitarias de presencia-ausencia solamente (matrices de incidencia).

A continuación ofrecemos un ejemplo de la función “sampleSpecies” con algunos datos reales de una comunidad de plantas de pradera en la República Checa. Originalmente, este conjunto de datos contiene 61 especies en 12 parcelas, y los datos de rasgos que acompañan a las especies. Los datos de la comunidad y los datos de los rasgos se almacenan en los elementos de un objeto “lista” llamado ohrazeni, por el nombre de la ubicación del prado. Estos datos forman parte del paquete traitor, por lo que primero tenemos que cargarlos con la función “data”. A efectos de demostración, reducimos las 60 especies originales a 10 especies.


ohrazeni$vegetationdata <- ohrazeni$vegetationdata[, 1:10]
ohrazeni$traitdata <- ohrazeni$traitdata[1:10, ]
ohrazeni
#> $vegetationdata
#>       sp1    sp2   sp3   sp4   sp5   sp6    sp7   sp8   sp9  sp10
#> 1   0.010  0.000 3.863 0.000 0.868 2.103  0.000 3.135 0.021 0.029
#> 2  27.675  0.016 0.665 0.232 1.387 0.068  5.928 0.469 0.012 0.000
#> 3   0.240  0.000 0.136 0.000 0.377 1.113  0.000 1.492 0.000 0.000
#> 4  77.943  0.000 0.129 0.000 0.305 0.862  0.071 0.249 0.000 0.000
#> 5   0.347 28.678 0.992 0.410 1.611 0.002  0.000 0.000 0.000 0.000
#> 6   0.019  1.245 1.639 0.000 0.537 1.130  0.000 0.727 0.000 0.000
#> 7  24.808 10.784 0.000 0.000 0.048 1.360  0.173 0.039 0.000 0.000
#> 8   0.000  3.035 5.603 0.000 0.130 3.326  0.000 3.573 0.055 0.000
#> 9   0.000 27.886 0.000 0.000 0.000 0.000 57.711 0.000 0.000 0.000
#> 10  2.185  0.000 0.080 0.000 0.184 1.185  1.035 0.838 0.000 0.000
#> 11  0.933  5.020 0.000 0.000 0.086 0.970  0.000 1.401 0.000 0.000
#> 12 11.282  0.000 0.000 0.000 1.985 0.359  0.225 0.000 0.000 0.000
#> 
#> $traitdata
#>      lifeform growthform leafposition   sla  ldmc height seedweight
#> sp1      hemi      grass      semiros 28.30 275.9   60.5        0.0
#> sp2      hemi      grass      semiros 27.70 294.3   39.5        0.1
#> sp3      hemi       forb        erosu 23.80 227.4  105.0        0.2
#> sp4      hemi       forb      semiros 30.40 166.2   45.0        1.3
#> sp5      hemi       forb      semiros 28.40 187.4   75.0        0.7
#> sp6      hemi      grass      semiros 23.30 322.0   27.5        0.5
#> sp7      hemi       forb      semiros 25.00 217.1   65.0        0.6
#> sp8      hemi      grass      semiros 18.90 254.5   32.5        0.1
#> sp9      hemi      sedge      semiros 17.93 365.5   25.0        0.7
#> sp10     hemi       forb      semiros 26.10 161.5   90.0        0.6

En este raro caso, no hay NAs en los datos de los rasgos. Para simular un conjunto de datos más realista, creamos algunos NA para uno de los rasgos sustituyendo los valores existentes por NA. Existen en realidad dos formas de especificar los datos que faltan, a través de NA o como se muestra arriba a través de un vector de 0 y 1, donde 0 significa que no hay datos de rasgos disponibles. En este caso utilizamos directamente el vector de valores de rasgos, en el que las especies con datos de rasgos ausentes se definen mediante los NA. Si en ambos conjuntos de datos (datos de rasgo y de comunidad) se dan nombres de especies, éstos se emparejan automáticamente mediante dichos nombres. Si no se nombran como tales, se supone que las especies de los datos de rasgos y de la comunidad están exactamente en el mismo orden. Si no lo están, sampleSpecies indicará la especie incorrecta a muestrear. Como ya habrá observado en los ejemplos anteriores, la función emitirá una advertencia recordándole esto, cuando no se den nombres de especies.


sla_with_NA <- ohrazeni$traitdata$sla
sla_with_NA[sample(1:10, 4)] <- NA
names(sla_with_NA) <- rownames(ohrazeni$traitdata)
sampleSpecies(ohrazeni$vegetationdata, sla_with_NA, thresh = 0.9,
              sequential = FALSE)
#> $Original
#> $Original$com_available
#>       sp2   sp4   sp5    sp7   sp8  sp10
#> 1   0.000 0.000 0.868  0.000 3.135 0.029
#> 2   0.016 0.232 1.387  5.928 0.469 0.000
#> 3   0.000 0.000 0.377  0.000 1.492 0.000
#> 4   0.000 0.000 0.305  0.071 0.249 0.000
#> 5  28.678 0.410 1.611  0.000 0.000 0.000
#> 6   1.245 0.000 0.537  0.000 0.727 0.000
#> 7  10.784 0.000 0.048  0.173 0.039 0.000
#> 8   3.035 0.000 0.130  0.000 3.573 0.000
#> 9  27.886 0.000 0.000 57.711 0.000 0.000
#> 10  0.000 0.000 0.184  1.035 0.838 0.000
#> 11  5.020 0.000 0.086  0.000 1.401 0.000
#> 12  0.000 0.000 1.985  0.225 0.000 0.000
#> 
#> $Original$available_pool
#> [1] 0.4861936
#> 
#> $Original$available_plot
#>  [1] 0.402034101 0.220344563 0.556581298 0.007855805 0.958146067 0.473664338
#>  [7] 0.296785983 0.428571429 1.000000000 0.373524605 0.773721760 0.159555267
#> 
#> 
#> $Scenario1
#> $Scenario1$sample_species
#> [1] "sp1"
#> 
#> $Scenario1$com_available
#>       sp1    sp2 sp3   sp4   sp5 sp6    sp7   sp8 sp9  sp10
#> 1   0.010  0.000  NA 0.000 0.868  NA  0.000 3.135  NA 0.029
#> 2  27.675  0.016  NA 0.232 1.387  NA  5.928 0.469  NA 0.000
#> 3   0.240  0.000  NA 0.000 0.377  NA  0.000 1.492  NA 0.000
#> 4  77.943  0.000  NA 0.000 0.305  NA  0.071 0.249  NA 0.000
#> 5   0.347 28.678  NA 0.410 1.611  NA  0.000 0.000  NA 0.000
#> 6   0.019  1.245  NA 0.000 0.537  NA  0.000 0.727  NA 0.000
#> 7  24.808 10.784  NA 0.000 0.048  NA  0.173 0.039  NA 0.000
#> 8   0.000  3.035  NA 0.000 0.130  NA  0.000 3.573  NA 0.000
#> 9   0.000 27.886  NA 0.000 0.000  NA 57.711 0.000  NA 0.000
#> 10  2.185  0.000  NA 0.000 0.184  NA  1.035 0.838  NA 0.000
#> 11  0.933  5.020  NA 0.000 0.086  NA  0.000 1.401  NA 0.000
#> 12 11.282  0.000  NA 0.000 1.985  NA  0.225 0.000  NA 0.000
#> 
#> $Scenario1$available_pool
#> [1] 0.9229118
#> 
#> $Scenario1$available_plot
#>  [1] 0.4030312 0.9795622 0.6280524 0.9875438 0.9689763 0.4772513 0.9634526
#>  [8] 0.4285714 1.0000000 0.7702924 0.8846611 0.9740813
#> 
#> 
#> $Scenario2
#> $Scenario2$sample_species
#> $Scenario2$sample_species[[1]]
#> [1] "sp3" "sp6"
#> 
#> $Scenario2$sample_species[[2]]
#> [1] "sp1"
#> 
#> $Scenario2$sample_species[[3]]
#> [1] "sp1" "sp6"
#> 
#> $Scenario2$sample_species[[4]]
#> [1] "sp1"
#> 
#> $Scenario2$sample_species[[5]]
#> [1] "none"
#> 
#> $Scenario2$sample_species[[6]]
#> [1] "sp3" "sp6"
#> 
#> $Scenario2$sample_species[[7]]
#> [1] "sp1"
#> 
#> $Scenario2$sample_species[[8]]
#> [1] "sp3" "sp6"
#> 
#> $Scenario2$sample_species[[9]]
#> [1] "none"
#> 
#> $Scenario2$sample_species[[10]]
#> [1] "sp1" "sp6"
#> 
#> $Scenario2$sample_species[[11]]
#> [1] "sp1" "sp6"
#> 
#> $Scenario2$sample_species[[12]]
#> [1] "sp1"
#> 
#> 
#> $Scenario2$com_available
#>       sp1    sp2   sp3   sp4   sp5   sp6    sp7   sp8 sp9  sp10
#> 1      NA  0.000 3.863 0.000 0.868 2.103  0.000 3.135  NA 0.029
#> 2  27.675  0.016    NA 0.232 1.387    NA  5.928 0.469  NA 0.000
#> 3   0.240  0.000    NA 0.000 0.377 1.113  0.000 1.492  NA 0.000
#> 4  77.943  0.000    NA 0.000 0.305    NA  0.071 0.249  NA 0.000
#> 5      NA 28.678    NA 0.410 1.611    NA  0.000 0.000  NA 0.000
#> 6      NA  1.245 1.639 0.000 0.537 1.130  0.000 0.727  NA 0.000
#> 7  24.808 10.784    NA 0.000 0.048    NA  0.173 0.039  NA 0.000
#> 8      NA  3.035 5.603 0.000 0.130 3.326  0.000 3.573  NA 0.000
#> 9      NA 27.886    NA 0.000 0.000    NA 57.711 0.000  NA 0.000
#> 10  2.185  0.000    NA 0.000 0.184 1.185  1.035 0.838  NA 0.000
#> 11  0.933  5.020    NA 0.000 0.086 0.970  0.000 1.401  NA 0.000
#> 12 11.282  0.000    NA 0.000 1.985    NA  0.225 0.000  NA 0.000
#> 
#> $Scenario2$available_pool
#> [1] 0.9846352
#> 
#> $Scenario2$available_plot
#>         1         2         3         4         5         6         7         8 
#> 0.9969090 0.9795622 0.9594997 0.9875438 0.9581461 0.9964131 0.9634526 0.9965017 
#>         9        10        11        12 
#> 1.0000000 0.9854730 1.0000000 0.9740813

Disimilitud de rasgos

En este ejercicio aprenderemos a calcular e interpretar la disimilitud de rasgos entre pares de organismos. La mayoría de los ejemplos se centrarán en la estimación de la disimilitud entre pares de especies, pero el mismo enfoque puede utilizarse para estimar la disimilitud entre individuos (como mostraremos en un singular ejemplo más adelante).

Datos, paquetes y funciones adicionales

Trabajaremos primero con algunos datos inventados. A continuación, utilizaremos algunos rasgos de estudiantes, obtenidos en uno de nuestros cursos anteriores (“students_traits.txt”). Por último, consideraremos algunos datos de campo (“spxtNESpain.txt”) con rasgos de codificación difusa de unas 130 especies recogidas a lo largo de un gradiente climático en el NE de España (de Bello et al., 2006).

# install.packages(c("FD", "NbClust", "gtools"))
library(FD)
library(NbClust)
library(gtools)

url.ch3 <- "./0Data/Trait_Based_Ecology/chapter3/"
  
students <- read.table(paste(url.ch3,"students_traits.txt", sep = ""),
                       header = TRUE)
spxtNESpain <- read.table(paste(url.ch3,"spxtNESpain.txt", sep = ""),
                          header = TRUE)
head(students)

head(spxtNESpain)


# Se requiere cargar la siguiente función:  
source(paste(url.ch3,"trova.r", sep = ""))

Cálculo de la disimilitud de rasgos

Para simplificar comenzamos por crear una matriz especie por rasgos a mano. Esta incluirá 3 rasgos: (1) un rasgo cuantitativo, tamaño corporal (t1); (2) un rasgo cualitativo transformado a binario, si las especies son carnívoras o no (t2; 1=sí, 0=no); (3) y un rasgo categórico - o factor - con múltiples niveles (el color de las especies (t3) utilizando codificación difusa. A continuación explicamos con más detalle el razonamiento que hay detrás de dicha codificación difusa.


bodysize <- c(10, 20, 30, 40, 50, NA, 70)
carnivory <- c(1, 1, 0, 1, 0,1, 0)
red <- c(1, 0, 0.5, 0, 0.2, 0, 1)
yellow <- c(0, 1, 0, 0, 0.3, 1, 0)
blue <- c(0, 0, 0.5,1, 0.5, 0, 0)
colors.fuzzy <- cbind(red, yellow, blue)
names(bodysize) <- paste("sp", 1:7, sep="")
names(carnivory) <- paste("sp", 1:7, sep="")
rownames(colors.fuzzy) <- paste("sp", 1:7, sep="")
tall <- as.data.frame(cbind(bodysize, carnivory, colors.fuzzy))
tall

Comencemos ahora por calcular la disimilitud a mano, sin utilizar ninguna función existente. Lo haremos para cada rasgo por separado, considerando primero el rasgo binario, por ejemplo, especies carnívoras frente a no carnívoras (1 y 0 respectivamente en el rasgo carnivoría en la Fig. 3.3). La disimilitud entre dos especies viene dada, en este caso, simplemente por si las dos especies comparten el mismo nivel de carnivoría o no. Por ejemplo, las especies 1 y 2 son ambas carnívoras, por lo que la disimilitud entre ellas es cero. Esto puede calcularse simplemente como la diferencia del valor del rasgo de la especie 1 (que es 1) y la especie 2 (también 1), como $abs(1-1)=0$ (nótese que expresamos esto como diferencia absoluta porque las disimilitudes no pueden ser negativas). Por el contrario las especies 1 y 3 son completamente diferentes en este rasgo, una es carnívora y la otra no. En ese caso $abs(1-0)=1$ lo que implica que las dos especies son 100% diferentes entre sí. Siguiendo esto, podemos calcular la disimilitud para todos los pares de especies utilizando sólo la carnivoría. Esto puede hacerse simplemente considerando la función “dist” (primer línea de código, debajo). Para algunas funciones necesitamos la disimilitud de rasgos en tales datos de forma “triangular”. Para otras (por ejemplo, “mpd” en “picante”) necesitamos transformar tal objeto “triangular” en una matriz (segunda línea de código).


dist(tall$carnivory)
#>   1 2 3 4 5 6
#> 2 0          
#> 3 1 1        
#> 4 0 0 1      
#> 5 1 1 0 1    
#> 6 0 0 1 0 1  
#> 7 1 1 0 1 0 1
as.matrix(dist(carnivory))
#>     sp1 sp2 sp3 sp4 sp5 sp6 sp7
#> sp1   0   0   1   0   1   0   1
#> sp2   0   0   1   0   1   0   1
#> sp3   1   1   0   1   0   1   0
#> sp4   0   0   1   0   1   0   1
#> sp5   1   1   0   1   0   1   0
#> sp6   0   0   1   0   1   0   1
#> sp7   1   1   0   1   0   1   0

Observen que al codificar el rasgo “carnivoría” entre 0 y 1, todas las disimilitudes entre las especies se restringen entre 0 y 1 (indicando que las especies son iguales entre sí o completamente diferentes). Es importante observar que la diagonal de matrices de disimilitud es, en la mayoría de los casos, igual a cero. Esto es así porque en la mayoría de los enfoques se considera que la disimilitud de una especie consigo misma es igual a cero, lo que significa que una especie es igual (es decir, no diferente) a sí misma. Lo ideal sería incluir la disimilitud entre individuos dentro de una especie, pero la mayoría de los algoritmos no consideran este caso.

Consideremos ahora un rasgo cuantitativo, el tamaño corporal. Intuitivamente, la disimilitud entre estudiantes de 1,80 m y de 1,70 m de altura es de 10 cm. Análogamente, la disimilitud del tamaño corporal entre especies 1 y 2 será simplemente la diferencia en sus valores de rasgo observados $abs(10-20)=10$; entre las especies 1 y 3 $abs(10-30)=20$, y así sucesivamente. Para comparar estos valores de disimilitud con los resultados de otros rasgos, por ejemplo sobre la carnivoría (que se expresa en relación con su valor máximo posible, 1), todos los rasgos deben ser escalados a unidades similares. A este fin, en el caso de los rasgos cuantitativos, hay que definir cuál es la máxima disimilitud posible (primer línea de código, debajo). Otros enfoques utilizan la función “scale”, que calcula la media (centro) de un rasgo cuantitativo y las diferencias con respecto a dicha media. Luego, los valores se estandarizan por la desviación estándar del rasgo. Esto sólo puede aplicarse a rasgos cuantitativos (segunda línea de código). Los dos enfoques descritos (estandarización entre 0 y 1 versus el escalamiento) son generalmente comparables. La ventaja del primero es que puede utilizarse más fácilmente para comparar rasgos cuantitativos y cualitativos (categóricos).


tall$bodysize / max(tall$bodysize, na.rm=T) # Noten que tenemos
#> [1] 0.1428571 0.2857143 0.4285714 0.5714286 0.7142857        NA 1.0000000
                                            # un NA a tener en cuenta
(tall$bodysize - mean(tall$bodysize, na.rm=T)) / sd(tall$bodysize,
                                                    na.rm = T) # a mano
#> [1] -1.2344268 -0.7715167 -0.3086067  0.1543033  0.6172134         NA  1.5430335
scale(tall$bodysize) # con la función integrada  
#>            [,1]
#> [1,] -1.2344268
#> [2,] -0.7715167
#> [3,] -0.3086067
#> [4,]  0.1543033
#> [5,]  0.6172134
#> [6,]         NA
#> [7,]  1.5430335
#> attr(,"scaled:center")
#> [1] 36.66667
#> attr(,"scaled:scale")
#> [1] 21.60247

Uno de los enfoques más utilizados para calcular la disimilitud de los rasgos, la distancia de Gower, sigue el primer enfoque. Considera, para los rasgos cuantitativos, la máxima disimilitud para un rasgo dado en un conjunto de datos para estandarizar los rasgos. En nuestro caso, la máxima disimilitud de tamaños corporales viene dada por la diferencia entre la especie más grande y la más pequeña, esto es, $abs(70-10)=60$. Una vez que conocemos la máxima disimilitud, podemos simplemente dividir la disimilitud entre las especies, en su escala original, por dicho valor de disimilitud máxima en el conjunto de datos.


dist(tall$bodysize) / max(dist(tall$bodysize), na.rm = T) # Gower a mano  
#>           1         2         3         4         5         6
#> 2 0.1666667                                                  
#> 3 0.3333333 0.1666667                                        
#> 4 0.5000000 0.3333333 0.1666667                              
#> 5 0.6666667 0.5000000 0.3333333 0.1666667                    
#> 6        NA        NA        NA        NA        NA          
#> 7 1.0000000 0.8333333 0.6666667 0.5000000 0.3333333        NA

Notar que tenemos una especie con un valor faltante, por lo que la disimilitud entre esta especie (sp6) y el resto de las especies no puede ser calculada al menos para este rasgo. Esto es un problema, ya que la mayoría de las funciones que utilizaremos para calcular la diversidad funcional no pueden manejar los NA. Veremos más adelante que esto no nos impedirá calcular la disimilitud basada en todos los rasgos juntos, si tenemos información suﬀciente para otros rasgos. Escalar la disimilitud entre 0 y 1 no sólo se recomienda para comparar los rasgos cuantitativos con otro tipo de rasgos, sino que facilita algunas interpretaciones de la disimilitud y hace más robusto el cálculo de la diversidad funcional $\alpha$ y $\beta$. La distancia de Gower, que toma valores entre 0 a 1 para los rasgos cuantitativos, representa lo relativamente diferentes que son dos especies en comparación con el valor máximo de disimilitud posible. Por ejemplo, la disimilitud entre las especies 1 y 3 = 0,33 en términos de tamaño corporal, implica que las especies 1 y 3 son un 33% diferentes de la disimilitud máxima posible (que es de 60 cm). La gran ventaja es que podemos comparar rasgos independientemente de la unidad (por ejemplo, metros frente a gramos) porque la disimilitud de cada rasgo estaría estandarizada a su valor máximo. Del mismo modo, podemos comparar rasgos cuantitativos con otro tipo de rasgos (entre las especies 1 y 2 es 0,17 para el tamaño corporal 1 para la carnivoría). Eventualmente, es posible combinar estas dos disimilitudes. Utilizando la disimilitud de Gower, esto se expresa simplemente como un promedio de las dos disimilitudes, utilizando la media aritmética ($(0.17+1)/2=0.58$) o geométrica (distancia Euclidiana) ($sqrt(0.17^{2+1/2}=1.01$). Para la disimilitud de Gower, el paquete FD incluye una función útil, gowdis, que calcula por defecto la media aritmética.


library(FD)
gowdis(tall[, c("bodysize", "carnivory")]) # combinando 2 rasgos
#>            sp1        sp2        sp3        sp4        sp5        sp6
#> sp2 0.08333333                                                       
#> sp3 0.66666667 0.58333333                                            
#> sp4 0.25000000 0.16666667 0.58333333                                 
#> sp5 0.83333333 0.75000000 0.16666667 0.58333333                      
#> sp6 0.00000000 0.00000000 1.00000000 0.00000000 1.00000000           
#> sp7 1.00000000 0.91666667 0.33333333 0.75000000 0.16666667 1.00000000

# A mano sería:
(dist(tall$carnivory)+dist(tall$bodysize)/max(dist(tall$bodysize), na.rm=T))/2
#>            1          2          3          4          5          6
#> 2 0.08333333                                                       
#> 3 0.66666667 0.58333333                                            
#> 4 0.25000000 0.16666667 0.58333333                                 
#> 5 0.83333333 0.75000000 0.16666667 0.58333333                      
#> 6         NA         NA         NA         NA         NA           
#> 7 1.00000000 0.91666667 0.33333333 0.75000000 0.16666667         NA

La única diferencia en estos dos enfoques es que la función gowdis no proporciona NAs. Esto se debe a que en lugar de calcular la disimilitud media entre los dos rasgos, para aquellas especies con NA, sólo se utiliza la disimilitud de los rasgos sin NAs (para sp6 sólo se utiliza la carnivoría en gowdis). Regresaremos a esto más adelante.

Algunos problemas con la distancia de Gower

Hay que tener en cuenta que la combinación de rasgos cuantitativos y cualitativos puede plantear varios problemas, que a menudo no se reconocen en la literatura. Imaginemos un caso, resumido en el siguiente script de R, en el que tenemos 20 especies y dos rasgos, un rasgo cuantitativo (en el que los valores entre 1 y 100 se extraen al azar) y un rasgo binario con valores 0 y 1 asignados aleatoriamente. Si calculamos la disimilitud para cada rasgo y luego combinamos estas dos disimilitudes con una media simple (o incluso con la distancia euclidiana, es decir, la media geométrica), la disimilitud resultante se verá mucho más afectada por el rasgo binario que por el cuantitativo (figura 1). En otras palabras, la contribución del rasgo binario a la disimilitud global será desproporcionada. Esto se debe a que los valores de disimilitud con un rasgo binario son más extremos, presentan una mayor varianza (figura 1, Pavoine et al., 2009). Las disimilitudes extremas son más comunes en el rasgo binario: en nuestro ejemplo, todas las especies que pertenecen a un nivel de carnivoría diferente tendrán una disimilitud de 1, mientras que para el tamaño corporal sólo un par de especies (la especie 1 y la especie 7) tienen este valor de disimilitud extrema. Como consecuencia, la disimilitud combinada estará mucho más correlacionada con la disimilitud calculada utilizando sólo el rasgo binario (Pearson $r>0,85$, en múltiples corridas), en comparación con la correlación con el rasgo cuantitativo solamente (Pearson $r<0,5$). Tengan en cuenta que, dado que el script de abajo se basa en valores aleatorios, estos resultados podrían cambiar cuando ejecuten el script).


quantitative <- sample(1:100, 20)
binary <- sample(0:1, 20, replace = T)
dissim.quant <- dist(quantitative) / max(dist(quantitative))
dissim.bin <- dist(binary)
# comb.dissim <- (dissim.quant+dissim.bin)/2 # distancia de gower a mano
comb.dissim <- gowdis(cbind(quantitative, binary)) # hace lo mismo con gowdis

cor(dissim.quant, comb.dissim)
#> [1] 0.3427852
cor(dissim.bin, comb.dissim)
#> [1] 0.9019734

par(mfrow = c(1, 3))
hist(dissim.quant, main = "quantitative", col = "grey")
hist(dissim.bin, main = "binary", col = "grey")
hist(comb.dissim, main = "combined", col = "grey")

$\label{fig:figs} Figura 1. Distribución de los valores de dos rasgos> cuantitativo (izquierda), y binario (centro). En el panel de la derecha se muestra la distribución conjunta de ambos rasgos$

Figura 1. Distribución de los valores de dos rasgos> cuantitativo (izquierda), y binario (centro). En el panel de la derecha se muestra la distribución conjunta de ambos rasgos

Una posible solución es jugar con el argumento “w” en la función gowdis. Este argumento nos permite dar diferentes pesos a los rasgos. Entonces, al considerar la media de las disimilitudes entre los rasgos algunos rasgos podrían tener más peso (esto es, una media ponderada). Así, si damos más peso al rasgo cuantitativo (por ejemplo 2 veces el peso dado a los rasgos binarios), obtendremos una importancia más uniforme de los rasgos en la disimilitud combinada. En otras palabras, al disminuir el peso del rasgo categórico evitamos que este rasgo tenga una contribución desproporcionada a la disimilitud multi-rasgo.


comb.dissim.weighted <- gowdis(cbind(quantitative, binary), w = c(2, 1))

cor(dissim.quant, comb.dissim.weighted)
#> [1] 0.6371553
cor(dissim.bin, comb.dissim.weighted)
#> [1] 0.70569

Tengan en cuenta que más adelante - después de leer las siguientes secciones sobre los problemas con la distancia de Gower y los rasgos codificados difusos - presentamos una nueva función de R (una función alternativa a gowdis, la función gawdis, con “a”). Esta nueva función está diseñada para resolver el problema del escalamiento o estandarización cuando se combinan rasgos de distinta naturaleza (cunatitativos, categóricos binarios, difusos).

Rasgos codificados categóricos y difusos (“fuzzy”)

Consideremos ahora el tercer rasgo que hemos mencionado anteriormente, el color de la especie. ¿Cómo podemos transformar un rasgo con múltiples niveles “desordenados” (nominal) en un rasgo cuantitativo? ¿Es el rojo más grande - o mayor - que el azul? Ordenar un rasgo como éste es una tarea diﬀcil (¿qué es mayos o menor?). Quizás en este caso podríamos utilizar la concentración de algún pigmento para definir el color, resultando en una sola columna para definir el rasgo. Pero en la mayoría de las ocasiones tenemos información de rasgos expresada en múltiples niveles en los que es diﬀcil determinar qué es menor o mayor. Por ejemplo, las especies de algas pueden ser categorizadas por su rasgo morfológico y aspecto en coriáceas, foliosas, filamentosas, costrosas, etc. Esto no puede tratarse sólo como un valor, ya que cada una de estas estrategias se considera completamente diferente de otra. La fuente de alimentación es otro ejemplo, ya que suele haber varios niveles (los italianos comen espaguetis, los franceses croissant, los checos dumplings, los españoles tortilla, etc.). Para este tipo de rasgos, podríamos tener dos soluciones. Una es codificar el rasgo como un factor. Imaginemos un caso con 10 especies, cada una de las cuales es una hierba, o una legumbre o un maleza (3 primeras líneas del siguiente código).

Esto funciona básicamente como en el caso de la disimilitud con un rasgo binario. Si las especies están en la misma categoría, entonces la disimilitud es 0 y viceversa, si las especies están en diferentes categorías, entonces la disimilitud es 1. Así que podríamos tratar esto como un rasgo de tipo binario, pero como tenemos tres niveles (hierba, maleza, leguminosa) necesitaríamos 3 columnas (última línea de código, con función “acm.disjonctif”). Esta función crea las llamadas variables “dummy”, donde cada nivel de un factor tiene su propia columna. Observe que cada especie tiene un valor 1 para una columna y 0 para las demás. Así que la suma de las 3 columnas, para una especie dada, es 1 (la suma de las filas será igual a 1). OJO, la función gowdis, no necesita que realicemos este paso a mano, por lo que le dimos la información del rasgo directamente en el objeto “df”.


cat.trait <- sample(c("grass", "legume", "forb"), 10, replace = T)
df <- as.data.frame(cat.trait)
gowdis(as.data.frame(df))
#>    1 2 3 4 5 6 7 8 9
#> 2  0                
#> 3  1 1              
#> 4  1 1 1            
#> 5  1 1 0 1          
#> 6  0 0 1 1 1        
#> 7  0 0 1 1 1 0      
#> 8  1 1 1 0 1 1 1    
#> 9  1 1 1 0 1 1 1 0  
#> 10 1 1 1 0 1 1 1 0 0

acm.disjonctif(df)

El problema ocurre cuando una especie pertenece simultáneamente a diferentes categorías, por ejemplo en el caso (imposible) de que una especie pueda ser a la vez una hierba y una legumbre. ¿Pero qué pasa si los italianos no son los únicos que comen espaguetis y los checos también comen espaguetis? Lo mismo puede ocurrir con los colores anteriores, y con muchos rasgos categóricos. ¿Qué hacer aquí? Lo ideal es utilizar las variables “dummy” con un enfoque de codificación difusa. Por ejemplo, si el rasgo color tiene tres niveles principales (rojo, azul, amarillo, es decir, los tres colores primarios), entonces el rasgo puede resumirse en la matriz especie x rasgo mediante tres columnas, una para cada color, llamadas variables “dummy”. En las columnas, podemos entonces añadir 1 para una especie si pertenece completamente a una de las categorías, de las tres en este caso, a las que pertenece la especie. Por ejemplo, la especie 1 es roja, por lo que sumaremos 1 a la columna que se refiere al color rojo y 0 a las demás columnas. Para cada fila, la suma debe ser igual a 1, porque añadir 1 a una columna reflejaría que la especie tiene el 100% de un color determinado (ya que la especie 1 es 100% roja).


colors.fuzzy <- cbind(red, yellow, blue)
colors.fuzzy
#>      red yellow blue
#> [1,] 1.0    0.0  0.0
#> [2,] 0.0    1.0  0.0
#> [3,] 0.5    0.0  0.5
#> [4,] 0.0    0.0  1.0
#> [5,] 0.2    0.3  0.5
#> [6,] 0.0    1.0  0.0
#> [7,] 1.0    0.0  0.0

rowSums(colors.fuzzy)
#> [1] 1 1 1 1 1 1 1

La utilización de variables “dummy”, combinada con la codificación difusa, permite tener en cuenta que algunas especies tienen un valor intermedio entre los distintos niveles de los rasgos (es decir, que tienen valores que “se sitúan” entre las distintas categorías predefinidas del rasgo). Por ejemplo, una especie de color naranja podría codificarse como mitad roja y mitad amarilla. Entonces podríamos añadir 0,5 a la columna roja y 0,5 a la amarilla (es decir, la especie es 50% roja y 50% amarilla, lo que produce el color naranja). Si la especie es violeta, entonces sería mitad azul y mitad roja, como la especie 3 del ejemplo anterior. Estos casos intermedios son bastante frecuentes. Este enfoque de codificación difusa también puede utilizarse para los rasgos binarios, por ejemplo, algunas especies de plantas son a veces anuales, a veces perennes (o en el caso de la carnivoría mencionado anteriormente, una especie podría ser a veces carnívora y a veces no, una especie de vegetariana a tiempo parcial). Esto puede resolverse utilizando un rasgo binario que vaya de 0 a 1, pero considerando también valores intermedios, como 0,5 para las especies que son igualmente de una forma u otra. Imaginemos ahora una especie hipotética que a veces es carnívora y a veces no, y supongamos que codificamos esta especie como 0,5 para el rasgo “carnívoro”. Entonces, la disimilitud de esta especie con las especies puramente carnívoras o no carnívoras (que tienen 0 o 1 en carnivoría), sería $abs(1 - 0,5)$ o $abs(0 - 0,5)$, o sea, 0,5. Este valor potencial de 0,5 implicaría que las especies son disímiles en un 50%.

Veamos ahora cómo debe calcularse la disimilitud cuando se utilizan variables dummy, con codificación difusa o no. En primer lugar, supongamos que el amarillo, el rojo y el azul son colores completamente diferentes, para simplificar. Esta es una suposición que seguimos implícitamente al utilizar variables dummy. Entonces la disimilitud entre las especies 1 y 2 debería ser igual a 1 (véase el objeto colors.fuzzy anterior), porque cada especie “pertenece” a diferentes categorías del rasgo (se incluyen como 1 en las columnas roja y amarilla respectivamente, es decir, son 100% rojas y 100% amarillas, respectivamente). Por otro lado, la distancia entre las especies 1 y 3 debería ser de 0,5 porque la especie 3 es parcialmente roja (50%). Del mismo modo, la disimilitud entre las especies 1 y 5 debería ser de 0,8, porque las especies 1 y 5 sólo se solapan en el 20% de los casos. Hay varias formas posibles de calcular este tipo de distancias en R, pero hay que hacerlo con mucho cuidado. Lo más importante es que la mayoría de los usuarios no son conscientes de que aplicar la función comúnmente utilizada gowdis a la especie x rasgo llamada colors.fuzzy no funcionará bien, ya que no devolverá estos valores de disimilitud esperados. El problema es que la función no “entenderá” que las tres columnas que se refieren al color de la especie reflejan en realidad la misma información del mismo rasgo. Así que considerará a cada columna como un rasgo diferente y calculará un promedio de la disimilitud entre todas las columnas de la matriz especie x rasgo (primer línea de código). La segunda línea de código soluciona el problema. Pueden ver que la máxima disimilitud es 1, y es precisamente entre especies que no comparten ningún valor $> 0$ en una columna determinada.


gowdis(colors.fuzzy) # no da lo matriz de disimilitud deseada
#>           1         2         3         4         5         6
#> 2 0.6666667                                                  
#> 3 0.3333333 0.6666667                                        
#> 4 0.6666667 0.6666667 0.3333333                              
#> 5 0.5333333 0.4666667 0.2000000 0.3333333                    
#> 6 0.6666667 0.0000000 0.6666667 0.6666667 0.4666667          
#> 7 0.0000000 0.6666667 0.3333333 0.6666667 0.5333333 0.6666667
gowdis(colors.fuzzy) / max(gowdis(colors.fuzzy)) # ésta es la correcta
#>     1   2   3   4   5   6
#> 2 1.0                    
#> 3 0.5 1.0                
#> 4 1.0 1.0 0.5            
#> 5 0.8 0.7 0.3 0.5        
#> 6 1.0 0.0 1.0 1.0 0.7    
#> 7 0.0 1.0 0.5 1.0 0.8 1.0

El principal problema aparece si queremos calcular la disimilitud para el objeto “tall”, que incluye los 3 rasgos creados anteriormente (tamaño del cuerpo, carnivoría y los 3 colores). Al combinar los 3 rasgos para calcular una disimilitud conjunta, la función gowdis “pensará” que hay 5 rasgos. Al hacer el promedio de las disimilitudes para el rasgo único, dará 1/5 del peso a cada columna cuando debería dar 1/3 a los primeros dos rasgod y un peso de 3/5 (en lugar de 1/3) al rasgo color. Una tentación sería utilizar el argumento “w” en la función gowdis, pero esto no resolvería el problema. Lo mejor, en estos casos, es calcular las disimilitudes para los rasgos por separado y luego promediar estas disimilitudes entre los diferentes rasgos. Consideren que gowdis no funciona con vectores, sino que necesita una matriz o data.frame, como objeto que contiene los rasgos:


tall

dissim.bodysize <- gowdis(as.matrix(bodysize))
dissim.carnivory <- gowdis(as.matrix(carnivory))
dissim.colour <- gowdis(colors.fuzzy) / max(gowdis(colors.fuzzy))
(dissim.bodysize + dissim.carnivory + dissim.colour) / 3
#>           sp1       sp2       sp3       sp4       sp5       sp6
#> sp2 0.3888889                                                  
#> sp3 0.6111111 0.7222222                                        
#> sp4 0.5000000 0.4444444 0.5555556                              
#> sp5 0.8222222 0.7333333 0.2111111 0.5555556                    
#> sp6        NA        NA        NA        NA        NA          
#> sp7 0.6666667 0.9444444 0.3888889 0.8333333 0.3777778        NA

Valores faltantes (“NAs”)

Como vimos anteriormente, un problema importante en el cálculo de la disimilitud de rasgos es el caso de los valores de rasgos faltantes. En el caso anterior de la matriz “tall”, para el rasgo cuantitativo no se puede calcular la disimilitud entre la especie 6 y las demás especies (por los “NA”). La buena noticia es que la disimilitud aún puede calcularse teniendo en cuenta otros rasgos cuando, y sólo cuando, hay al menos un valor de rasgo disponible para la especie 6 (y la especie con la que se compara la especie 6). La disimilitud resultante que combina rasgos para la especie 6 es simplemente la media de la disimilitud para los rasgos sin NA. En el caso del ejemplo anterior, esto correspondería a la media entre la disimilitud de la carnivoría y el color solamente, excluyendo el tamaño del cuerpo (debido al NA). Entonces, podemos seguir utilizando la disimilitud resultante de otros rasgos y descartar, para las comparaciones con la especie 6, la información que falta sobre el tamaño del cuerpo. Por ejemplo, la distancia entre las especies 1 y 6 considerando el color y la carnivoría es 0,5, ya que la disimilitud basada en la carnivoría es igual a 0 y la disimilitud basada en el color es igual a 1 (vean “dissim.carnivory” y “dissim.colour”). Más adelante veremos cómo reemplazar los valores NA.

Grupos funcionales con rasgos de los estudiantes

Trabajemos ahora con los rasgos recogidos por los estudiantes. Algunes estudiantes asistieron a uno de nuestros cursos, en el que les pedimos que proporcionaran rasgos con el objetivo de dividirlos en diferentes “tipos”. Proporcionaron varios rasgos y aquí sólo utilizamos 3 de ellos. En este caso, les estudiantes tenían que decir cuánto alcohol consumen (en una escala del 1 al 10), cuántos kilómetros en bicicleta hacen a la semana y cuántas horas han tenido que viajar para venir al curso. En primer lugar, podemos ver que el último “rasgo” (por supuesto, no se trata de un rasgo funcional) no está distribuido de manera muy uniforme, porque la mayoría viven muy cerca de la universidad donde se imparte el curso. Así que quienes vienen de muy lejos parecerán muy diferentes. Sin embargo, si transformamos los datos en logaritmos obtendremos una distribución más “normal”, uniforme, y las disimilitudes estarán menos sesgadas (figura 2). Fíjate en que, cuando hacemos la transformación logarítmica, tenemos que añadir “+1”, ya que hay valores cero en el rasgo.


students

par(mfrow = c(2, 2))
hist(students$travel_dist, col = "grey", main = "trait values")
hist(log(students$travel_dist + 1), col = "grey", main = "trait values log")
hist(dist(students$travel_dist), col = "grey", main = "dissimilarities")
hist(dist(log(students$travel_dist + 1)), col = "grey",
     main = "dissimilarities log")

$\label{fig:figs} Figura 2. Distribución del rasgo travel_dist de los estudiantes: sin transformar (arriba e izquierda), transformación logarítmica (arriba, derecha). Los paneles de la fila inferior son análogos, pero calculados como distancias o disimilitudes$

Figura 2. Distribución del rasgo travel_dist de los estudiantes: sin transformar (arriba e izquierda), transformación logarítmica (arriba, derecha). Los paneles de la fila inferior son análogos, pero calculados como distancias o disimilitudes

El rasgo travel_dist no se distribuye normalmente, por lo que utilizamos la transformación logarítmica sustituyendo los valores originales en la matriz (figura 2). Aunque en algunos casos la normalidad no mejoró, al menos los valores no parecían tan dispersos después de la transformación. También podemos observar cómo se correlacionan los rasgos, para ver hasta qué punto los rasgos individuales podrían ser redundantes. En este caso, las correlaciones son relativamente bajas. Ahora podemos calcular la disimilitud de Gower entre los estudiantes, que en este caso es simple, ya que no hay rasgos codificados difusamente.


students$travel_dist <- log(students$travel_dist + 1)
cor(students)
#>                 alc_cons km_bike_week travel_dist
#> alc_cons      1.00000000    0.4029536 -0.07907633
#> km_bike_week  0.40295355    1.0000000  0.38094912
#> travel_dist  -0.07907633    0.3809491  1.00000000
distances.students <- gowdis(students)
# distances.students

El objeto resultante es relativamente grande y diﬀcil de visualizar. ¿Cómo podríamos visualizar cuánto se parecen/diferencian los estudiantes entre sí? Podemos utilizar el PCoA para ello, ya que se basa en una matriz de disimilitud (figura 3). Obsérvese que en este caso podríamos haber utilizado incluso PCA, porque todos los rasgos eran cuantitativos y no hay valores “NA”, pero a menudo no se cumplen estas condiciones. Así que vamos a ejecutar el PCoA. A veces, como en este caso, utilizar la raíz cuadrada de la distancia mejora las características de la disimilitud y la función “dudi.pco” es más amigable (no da mensajes de advertencia).


pcoa.all <- dudi.pco(sqrt(distances.students), scannf = F, nf = 4)
scatter(pcoa.all)

$\label{fig:figs} Figura 3. Ordenamiento de las distancias de los rasgos de les estudiantes con un análisis de correspondencia principal (PCoA)$

Figura 3. Ordenamiento de las distancias de los rasgos de les estudiantes con un análisis de correspondencia principal (PCoA)

sum(pcoa.all$eig[1:2]) / sum(pcoa.all$eig)
#> [1] 0.6402358
cor(pcoa.all$li[, 1:2], students)
#>      alc_cons km_bike_week travel_dist
#> A1 -0.1577509   -0.6710870  -0.9330768
#> A2 -0.9287152   -0.4986334   0.3018606

Los alumnos que están más cerca en la figura 3 tienen más rasgos similares (¡por ejemplo Ana y Joana deben ser muy parecidas entre ellas desde el punto de vista funcional!) Los dos primeros ejes del PCoA son los más importantes (sus autovalores - “eigenvalues” son más altos). La variabilidad capturada por los dos primeros ejes es $\approx 65$%. ¿Están estos ejes más relacionados con algunos rasgos que con otros? Podemos ver si la posición de los estudiantes a lo largo de un eje determinado (que está en el objeto “pcoa.all$li”) está correlacionada con algunos rasgos. Calculemos la correlación entre la posición de los alumnos en los dos primeros ejes del PCoA y los 3 rasgos. El primer eje, el que representa la mayor parte de las diferencias entre los estudiantes, representa en general su lugar de origen (las personas que viajan poco están en la parte izquierda/superior de la figura 3) y en parte cuánto utilizan la bicicleta (los que utilizan más la bicicleta están a la derecha). El segundo eje refleja el consumo de alcohol, la gente que está en la parte superior del gráfico bebe menos y también la distancia de desplazamiento (los que viajan más están en la parte inferior del gráfico) (figura 3). Como sólo tenemos rasgos cuantitativos, sin NA, también podríamos visualizar los resultados con un PCA que proporciona, más o menos, unos resultados comparables (figura 4).

par(mfrow=c(1, 2))
scatter(dudi.pco(sqrt(distances.students), scannf = F, nf = 4))
scatter(dudi.pca(students, scannf = F, nf = 4))

$\label{fig:figs} Figura 4. Comparación del ordenamiento de las distancias de los rasgos de les estudiantes con un análisis de correspondencia principal (PCoA)(izquierda) y un análisis de componentes principales (PCA)$

Figura 4. Comparación del ordenamiento de las distancias de los rasgos de les estudiantes con un análisis de correspondencia principal (PCoA)(izquierda) y un análisis de componentes principales (PCA)

Ahora podemos intentar ver si hay algunos “tipos” diferentes de estudiantes, como intentar hacer “grupos funcionales” de estudiantes basados en los rasgos considerados. Esto es fácil de calcular en R, pero los resultados serán generalmente bastante variables, porque hay muchas formas posibles de hacerlo. Podríamos utilizar una de las posibles opciones de agrupación, por ejemplo. Entonces podemos probar una de las muchas formas posibles de encontrar formalmente el número óptimo de grupos en dicho cluster (que maximice la disimilitud entre grupos y minimice la disimilitud dentro de los grupos). Una buena opción es la función “NbClust”, que ofrece muchos enfoques posibles. Esto nos dice que, si probamos cuántos grupos (entre 2 y 6) maximizarían la diferencia entre los estudiantes, ¡la respuesta es 5! Podemos visualizar qué estudiantes están en cada grupo explorando el objeto $Mejor.partición y también visualmente con un gráfico (figura 5).


par(mfrow=c(1, 2))
cluster.students <- hclust(distances.students, method = "ward.D2")
plot(cluster.students)
res.groups <- NbClust(diss = distances.students,
                      distance = NULL,
                      min.nc = 2, max.nc = 6,
                      method = "ward.D",
                      index = "silhouette")
#> 
#>  Only frey, mcclain, cindex, sihouette and dunn can be computed. To compute the other indices, data matrix is needed
res.groups
#> $All.index
#>      2      3      4      5      6 
#> 0.4309 0.5043 0.5105 0.5291 0.5282 
#> 
#> $Best.nc
#> Number_clusters     Value_Index 
#>          5.0000          0.5291 
#> 
#> $Best.partition
#>      Ana  Borbala Catarina  Eduardo    Elise Fernanda  Giorgia    Joana 
#>        1        2        1        3        2        3        2        1 
#> Johannes    Julia    Laura   Hannah  Martina  Norbert     Nuno    Sanne 
#>        4        2        5        2        1        4        3        5 
#> Veronika   Villem 
#>        2        2
plot(cluster.students)
rect.hclust(cluster.students, k = res.groups$Best.nc[1],
            border = 1:res.groups$Best.nc[1])

$\label{fig:figs} Figura 5. Dendrogramas del agrupamiento de les estudiantes por las distancias de sus rasgos (el dendrograma de la derecha muestra cinco grupos conformados)$

Figura 5. Dendrogramas del agrupamiento de les estudiantes por las distancias de sus rasgos (el dendrograma de la derecha muestra cinco grupos conformados)

También podemos visualizar los grupos utilizando “s.class”, con respecto al PCoA (figura 6):


par(mfrow = c(1, 2))
scatter(pcoa.all)
s.class(pcoa.all$li, as.factor(res.groups$Best.partition), cpoint = 1)

$\label{fig:figs} Figura 6. Conformación de los grupos en el plano de ordenamiento del PCoA)$

Figura 6. Conformación de los grupos en el plano de ordenamiento del PCoA)

La función “gawdis”

Para más información, consulte la página de CRAN del paquete gawdis

En esta guía paso a paso veremos cómo funciona la función “gawdis” y la aplicaremos utilizando datos básicos (disponibles en el paquete “FD”, o creados a continuación). La función “gawdis” es una extensión de la función clásica gowdis del paquete FD. Al final de la sección 3.3 vimos que la función gawdis (con “a”) fue diseñada como una alternativa a “gowdis”(con “o”) para resolver el problema del escalamiento o estandarización cuando se combinan rasgos de distinta naturaleza (cunatitativos, categóricos binarios, difusos). La función gawdis calcula la disimilitud entre unidades, normalmente especies, basándose en múltiples tipos de variables (cuantitativas o categóricas), normalmente rasgos de las especies. Por lo tanto, puede aplicarse para calcular la disimilitud de múltiples rasgos entre especies en estudios de ecología de rasgos funcionales, así como en otras aplicaciones. La disimilitud obtenida puede calcularse para obtener una contribución quasi-déntica de las variables individuales (por ejemplo, los rasgos) o del grupo de variables asociadas (en variables que reflejan información similar, como diversos rasgos de las hojas de una planta).

La disimilitud se calcula mediante un enfoque analítico o mediante iteraciones. La función toma prestados varios argumentos de “gowdis”, con otros adicionales que se describen a continuación. Incluye la opción de considerar variables difusas y “dummy” (esto es, múltiples columnas que caracterizan un único rasgo). Veamos los datos dummy$trait, que incluyen datos inventados para algunas especies y sus rasgos en el paquete FD. Los datos dummy.trait incluyen información de rasgos para 8 especies. Algunos rasgos son cuantitativos numéricos: “num1” y “num2”; otros son categóricos: los factores “fac1” y “fac2”; otros son rasgos semicuantitativos del tipo ordinal: “ord1” y “ord2”; y otros rasgos son binarios: bin1 y bin2. Para complicar las cosas, agunos rasgos tienen valores “NA”.

Cuando trabajamos con rasgos de diferente naturaleza (cuantitativos, binarios, etc) observamos que algunos rasgos, como los binarios, contribuyen más fuertemente a la disimilitud global (disimilitud multi-rasgo) que rasgos cuantitativos. Esto no es deseable, y exploramos la utilidad de pesar la importancia relativa de los rasgos con el argumento “w” en gowdis. En lugar de hacer una media simple, el argumento “w” nos posibilita hacer una media ponderada para emparejar la contribución de algunos rasgos. Podemos reducir el peso de rasgos binarios con el fin de reducir su contribución a la disimilitud multi-rasgo, y aumentar así la de los rasgos cuantitativos. Pero… ¿cómo sabemos qué valores tenemos que seleccionar para el argumento “w”, para que todos los rasgos tengan un peso similar? si sólo tenemos 2 rasgos quizá podamos probar varios valores de “w” y esperar encontrar una combinación decente. Pero la cuestión se complica si tenemos muchos rasgos. Aquí es donde la nueva función gawdis es útil. La función busca los mejores valores para el argumento “w”, para obtener una contribución equitativa de los rasgos individuales.

Para comenzar con un ejemplo simplificado tomamos sólo dos rasgos (“num2” y “bin2”) sin “NA” y vemos que gawdis, como gowdis, calcula una disimilitud para cada rasgo escalando las diferencias de valores entre cada par de especies a un valor entre 0 y 1. A diferencia de gowdis, gawdis calcula directamente una disimilitud multirasgo en la que la contribución de cada rasgo (proporcionada en la salida de la función, es decir, en “correls”) es básicamente idéntica. Esto es posible gracias a la estimación de un valor “w” ideal para cada rasgo (proporcionado en la salida de la función, “weights”).


library(FD)
# install.packages("gawdis")
library(gawdis)
dummy$trait

equalcont <- gawdis(dummy$trait[, c("num2", "bin2")])
#> [1] "Running w.type=analytic"
equalcont
#>            sp1        sp2        sp3        sp4        sp5        sp6
#> sp2 0.13584757                                                       
#> sp3 0.19924310 0.33509067                                            
#> sp4 0.42038371 0.39321419 0.61962681                                 
#> sp5 0.63773982 0.77358739 0.43849672 0.38037319                      
#> sp6 0.36226018 0.22641261 0.56150328 0.61962681 1.00000000           
#> sp7 0.55623127 0.69207884 0.35698817 0.29886465 0.08150854 0.91849146
#> sp8 0.18113009 0.04528252 0.38037319 0.43849672 0.81886991 0.18113009
#>            sp7
#> sp2           
#> sp3           
#> sp4           
#> sp5           
#> sp6           
#> sp7           
#> sp8 0.73736137
attr(equalcont, "correls") # contribución de cada rasgo
#>      num2      bin2 
#> 0.7377916 0.7377916
attr(equalcont, "weights") # pesos asignados para que contribuyan igual
#>      num2      bin2 
#> 0.6611248 0.3388752
# Para entender qué hace la función attr(), hagan el str() de equalcont,
# verán que el objeto equalcont, que es una matriz de distancias, tiene mucha
# informcación escondida muy útil.

Por supuesto, la función gawdis también puede utilizarse exactamente de la misma manera que la gowdis. Por ejemplo, podemos especificar que queremos tener un peso similar (y por tanto diferente contribución) para los diferentes rasgos. Este es el valor por defecto en gowdis y en gawdis esto se establece ya sea usando w.type =“equal” o por w.type =“user” y luego diciendo que “W” es el mismo para todos los rasgos (vean la correlaciones que indican idéntiicos comportamientos entre gowdis y gawdis en la figura 7).


ex1 <- gowdis(dummy$trait)
ex1.gaw1 <- gawdis(dummy$trait, w.type = "equal")
#> [1] "Running w.type=equal"
ex1.gaw2 <- gawdis(dummy$trait, w.type = "user", W = rep(1, 8))
#> [1] "Running w.type=user"

par(mfrow=c(1, 2))
plot(ex1, ex1.gaw1)
abline(0, 1)
plot(ex1, ex1.gaw2)
abline(0, 1)

$\label{fig:figs} Figura 7. Correlaciones para mostrar resultados idénticos de la función gowdis y gawdis variando el argumento w.type$

Figura 7. Correlaciones para mostrar resultados idénticos de la función gowdis y gawdis variando el argumento w.type

Las soluciones proporcionadas por la función se pueden obtener de dos maneras. La primera es mediante una solución analítica, utilizando fórmulas puramente matemáticas (véase el documento principal y el material suplementario correspondiente). Este es el enfoque que se utiliza por defecto en la función, tal y como hemos utilizado para crear el objeto equalcont anteriormente. En caso de duda, este enfoque se puede establecer utilizando w.type=“analytic”. Tengan en cuenta que, lamentablemente, esta solución no se puede utilizar cuando hay valores perdidos (NA).

Cuando hay NAs, incluso si se establece el argumento w.type=“analytic” la función aplicará un segundo enfoque, basado en iteraciones, es decir, seguirá probando varias soluciones hasta encontrar una adecuada (se basa en realidad en un enfoque de optimización de la aptitud, que progresivamente mejora la salida de los resultados). Esto llevará algún tiempo, ya que se prueban múltiples alternativas de forma secuencial. Este enfoque se implementa utilizando w.type =“optimized”. El tiempo de ejecución dependerá, sobre todo, del número de especies y del número de iteraciones, que se establece mediante el argumento “opti.maxiter”, fijado por defecto en 300. A continuación aplicamos ambos enfoques a un subconjunto de datos de dummy$trait, sin NAs, y verificamos mediante una correlación la similitud de ambos enfoques (figura 8).


analytical <- gawdis(dummy$trait[, c(2, 4, 6, 8)], w.type = "analytic")
#> [1] "Running w.type=analytic"
# No es necesario añadir el argumento w.type, es el argumento
# utilizado por defecto, si no se especifica.  

attr(analytical, "correls")
#>      num2      fac2      ord2      bin2 
#> 0.5350139 0.5350139 0.5350139 0.5350139
attr(analytical, "weights")
#>      num2      fac2      ord2      bin2 
#> 0.2422698 0.2467646 0.3575842 0.1533815

iterations <- gawdis(dummy$trait[, c(2, 4, 6, 8)],
                     w.type = "optimized",
                     opti.maxiter = 100)
#> [1] "Running w.type=optimized"
# aquí hemos utilizado "sólo" 100 iteraciones, para acelerar el proceso
# y porque con tan pocas especies esto es probablemente más que suficiente
attr(iterations, "correls")
#>      num2      fac2      ord2      bin2 
#> 0.5352085 0.5347223 0.5343517 0.5361612
attr(iterations, "weights")
#>      num2      fac2      ord2      bin2 
#> 0.2424975 0.2465992 0.3569872 0.1539161
plot(analytical, iterations)

$\label{fig:figs} Figura 8. Correlaciones para comparar la similitud de los resultados de los enfoques analítico y de optimización$

Figura 8. Correlaciones para comparar la similitud de los resultados de los enfoques analítico y de optimización

Se advierten aquí algunas ligeras diferencias en los resultados, por ejemplo en lo que respecta a las correlaciones entre la disimilitud de los rasgos individuales y el multirasgo (es decir, las “correlaciones”) y las ponderaciones finalmente utilizadas para cada rasgo (“pesos”)(figura 8). Los resultados finales, los valores de disimilitud mostrados en el gráfico, también varían un poco. Esto se debe a que las iteraciones se basan simplemente en un enfoque de prueba-error, que mejora con más iteraciones, pero que puede “perder” la solución matemática perfecta, y sólo acercarse a ella. Por lo demás, podemos estar seguros de que los resultados obtenidos con las iteraciones se acercan bastante a los resultados ideales.

Agrupación de rasgos con “gawdis”

En algunos casos, los conjuntos de datos que estamos considerando tienen múltiples variables que reflejan una información parcialmente solapada o redundante. Este es el caso del conjunto de datos “tussock.trait”, del paquete “FD”. Hay cinco rasgos foliares, probablemente medidos en las mismas hojas, y bastante correlacionados entre sí. Si muchos rasgos están muy correlacionados, proporcionan una información similar y tendrán una influencia prominente en la disimilitud de múltiples rasgos. Para estos casos es recomendable crear grupos de rasgos (en un rato utilizando el argumento “groups” para ello). Primero reduciremos ligeramente el número de rasgos, por simplicidad, y porque el conjunto de datos incluye algunos rasgos categóricos muy desequilibrados (con pocas entradas para un número de categorías). Luego necesitamos transformar logarítmicamente algunos de los rasgos cuantitativos (la altura, la masa de semillas y hojas son los que necesitaron una transformación logarítmica en función de análisis previos). Consideren que si no aplicamos la transformación logarítmica los rasgos con una distribución anormal tendrán una mayor contribución por tener una mayor varianza.


dim(tussock$trait)
#> [1] 53 16
head(tussock$trait)

head(tussock$trait[, 3:7])

cor(tussock$trait[, 3:7], use = "complete")
#>             LDMC      leafN      leafP      leafS        SLA
#> LDMC   1.0000000 -0.6848213 -0.6682552 -0.5366135 -0.6417381
#> leafN -0.6848213  1.0000000  0.6052319  0.7546005  0.6384765
#> leafP -0.6682552  0.6052319  1.0000000  0.7703387  0.6212733
#> leafS -0.5366135  0.7546005  0.7703387  1.0000000  0.5966755
#> SLA   -0.6417381  0.6384765  0.6212733  0.5966755  1.0000000
tussock.trait <- tussock$trait[, c("height", "LDMC", "leafN","leafS",
                                   "leafP", "SLA", "seedmass", "raunkiaer",
                                   "pollination", "clonality", "growthform")]
tussock.trait.log <- tussock.trait # algunos rasgos necesitan una
                                   # transformación logarítmica  

# Una nueva matriz para almacenar los nuevos datos:
tussock.trait.log$height <- log(tussock.trait$height)
tussock.trait.log$seedmass <- log(tussock.trait$seedmass)
tussock.trait.log$leafS <- log(tussock.trait$leafS)
colnames(tussock.trait.log)
#>  [1] "height"      "LDMC"        "leafN"       "leafS"       "leafP"      
#>  [6] "SLA"         "seedmass"    "raunkiaer"   "pollination" "clonality"  
#> [11] "growthform"
tussock.trait.log

Ahora comparamos la distancia de Gower con “gowdis” y con “gawdis” (ésta última proporciona más opciones que “gowdis”, comparamos ambas salidas correlacionándolas) (figura 9). Calculamos también la contribución a la disimilitud para los rasgos foliares (que son los que estamos evaluando agrupar), para eso calculammos sus correlaciones con la disimilitud multi-rasgo (especificada en el argumento “correls”). Al final del “chunk” calculamos, con la ayuda del paquete vegan (haciendo una correlación entre dos matrices de disimilitud con la prueba de Mantel), la contribución a la dismilitud multi-rasgo del conjunto de rasgos foliares (o sea, del agrupamiento de esos rasgos foliares).



straightgowdis <- gowdis(tussock.trait.log) # la antigua función
straightgowdis.2 <- gawdis(tussock.trait.log, w.type = "equal", silent = T)
plot(straightgowdis, straightgowdis.2) # vemos que dan lo mismo.

$\label{fig:figs} Figura 9. Correlaciones para comparar la similitud entre los resultados de las funciones gowdis y gawdis aplicadas en un ejemplo real$

Figura 9. Correlaciones para comparar la similitud entre los resultados de las funciones gowdis y gawdis aplicadas en un ejemplo real


# Noten que lo anterior lo hice con fines comparativos, simplemente para
# mostrar que en este caso gawdis y gowdis dan lo mismo. La idea que en
# adelante usen gawdis (quees la función más moderna, y que permite pesar
# fácilmente y adecuadamente distintos tipos de variables). Lo que sigue
# rasgos (es para evaluar posibles agrupamientos de que es el objetivo de
# la sección).

cors.gow <- attr(straightgowdis.2, "correls") # contribución de cada rasgo

# Noten que los rasgos 2 a 6 son los rasgos foliares que estoy considerando
# agrupar. El "[12]" indica que agregaré en elemento en el lugar 12 al objeto
# cors.gow (verifiquen que antes de esto cors.gow tenía 11 elementos (los 11
# rasgos analizados). En otras palabras, estoy agregando la contribución del
# grupo foliar a cors.gow:

straightgowdis.3 <- gawdis(tussock.trait.log[, 2:6], w.type = "equal",
                           silent = T)

cors.gow[12] <- vegan::mantel(straightgowdis.2, straightgowdis.3,
                              na.rm = T)$statistic

# La expresión 'paquete::función' devuelve el valor de aplicar la función del
# paquete definido. En este caso, la función "mantel", del paquete "vegan",
# devuelve el estadístico Mantel, que resulta de correlacionar dos matrices
# de disimilitud (straightgowdis.2 y straightgowdis.3).  

# El código debajo simplemente da un nombre a la contribución del
# agrupamiento:
names(cors.gow)[12] <- "HOJAS"
cors.gow
#>      height        LDMC       leafN       leafS       leafP         SLA 
#>   0.2023004   0.5083692   0.4522899   0.4944064   0.4292520   0.3729605 
#>    seedmass   raunkiaer pollination   clonality  growthform       HOJAS 
#>   0.2604944   0.5941643   0.4528619   0.3237053   0.4668307   0.5889437

La mayor contribución (reflejada en los valores de “cors.gow”) es la del rasgo categórico “raunkiaer” (la forma de vida raunkiaer) y la combinación de rasgos de hoja (“HOJAS”). Esto se debe simplemente a que con w.type = “equal” la disimilitud multi-rasgo es un promedio de la disimilitud para rasgos individuales, por lo que los rasgos de hoja están representados 5/11 veces en este promedio, un efecto desproporcionado con respecto a otros rasgos. Podríamos decir que (aunque hay diferencias entre los rasgos de la hoja) constituyen el mismo TIPO de rasgo, pero contado 5 veces. Por eso es recomendable agrupar esos rasgos en el nuevo rasgo HOJAS. Esto evidencia que el problema no se resuelve utilizando un PCoA u otro análisis para reducción de la dimensionalidad o, en otras palabras, para reducir el número de rasgos y sintetizar los rasgos correlacionados en algún eje multivariante (como erróneamente es sugerido en mucha de la literatura actual). Una mejor solución puede obtenerse considerando que todos los rasgos de las hojas pertenecen a un grupo determinado. Así, gawdis calculará primero una disimilitud para todos los rasgos de la hoja juntos y luego intentará obtener para esta disimilitud combinada de la hoja el mismo peso que para otros rasgos. Veámoslo:


colnames(tussock.trait.log)
#>  [1] "height"      "LDMC"        "leafN"       "leafS"       "leafP"      
#>  [6] "SLA"         "seedmass"    "raunkiaer"   "pollination" "clonality"  
#> [11] "growthform"

# Como hay NAs en los datos, el procedimiento iterativo es el único posible.
gaw.groups <- gawdis(tussock.trait.log, w.type = "optimized",
                     opti.maxiter = 300,
                     groups.weight = T,
                     groups = c(1, 2, 2, 2, 2, 2, 3, 4, 5, 6, 7))
#> [1] "Running w.type=optimized on groups=c(1,2,2,2,2,2,3,4,5,6,7)"
#> [1] "Traits inside the group were weighted - optimized."

cors.gaw.gr <- attr(gaw.groups, "correls")
cors.gaw.gr[12] <- attr(gaw.groups, "group.correls")[2]
names(cors.gaw.gr)[12] <- "HOJAS"
cors.gaw.gr
#>      height        LDMC       leafN       leafS       leafP         SLA 
#>   0.4362648   0.3913038   0.3431186   0.3251978   0.2771238   0.2842367 
#>    seedmass   raunkiaer pollination   clonality  growthform       HOJAS 
#>   0.4297707   0.4439558   0.4483326   0.4496296   0.4483830   0.4435296

Si miramos el “cors.gaw.gr” podemos ver que individualmente los rasgos de las hojas tienen una contribución menor que otros rasgos. Sin embargo en combinación, como “HOJAS”, tienen una contribución comparable. Por supuesto, será difícil decidir cuándo y en qué caso debe definirse el grupo de rasgos. Pero creemos que si los rasgos se miden en el mismo órgano y proporcionan alguna información que se solapa, entonces, deberían considerarse como un grupo.

Codificación difusa y variables “dummy” con “gawdis”

Introducimos previamente que la función gawdis es particularmente útil en el caso de rasgos codificados como variables difusas o, como caso específico de esto, como variables “dummy”. Creemos primero una nueva matriz de rasgos, “tall”.


bodysize <- c(10, 20, 30, 40, 50, NA, 70)
carnivory <- c(1, 1, 0, 1, 0, 1, 0)
red <- c(1, 0, 0.5, 0, 0.2, 0, 1)
yellow <- c(0, 1, 0, 0, 0.3, 1, 0)
blue <- c(0, 0, 0.5, 1, 0.5, 0, 0)
colors.fuzzy <- cbind(red, yellow, blue)
names(bodysize) <- paste("sp", 1:7, sep = "")
names(carnivory) <- paste("sp", 1:7, sep = "")
rownames(colors.fuzzy) <- paste("sp", 1:7, sep = "")
tall <- as.data.frame(cbind(bodysize, carnivory, colors.fuzzy))
tall

El objeto “tall” incluye 3 rasgos para 7 especies, el tamaño del cuerpo (cuantitativo), la carnivoría (binaria/categórica) y el color. Este último está representado por 3 columnas, una para cada color básico (amarillo, rojo, azul). Por ejemplo, la primera especie (sp1), es roja, por lo que obtiene 1 en la columna “rojo” y 0 en las demás. Este tipo de variables, definidas por múltiples columnas, se llama variable dummy, y en este caso concreto de codificación difusa, porque el valor de cada columna puede ser diferente de 0 y 1, siendo la suma sobre las 3 columnas igual a 1. Para más información, haga clic aquí. Ciertamente no podemos aplicar la función “gowdis” a este tipo de matrices, porque la función “pensará” que hay un total de 5 rasgos, ya que la matriz contiene 5 columnas, y tratará cada una de las 3 columnas para los colores como rasgos independientes, resultando en una mayor contribución de este único rasgo. Además, la función se vuelve un poco “loca” con este tipo de variables.
Si queremos combinar este rasgo (el color, definido por 3 columnas en “tall”), necesitamos primero calcular la disimilitud interna para el color, y luego combinarla con los otros rasgos (por ejemplo con una media). Podemos resolver este problema con la función gawdis de forma sencilla. La solución se obtiene definiendo todas las columnas que pertenecen a un determinado rasgo a un determinado grupo, como vimos anteriormente, y luego estableciendo el argumento fuzzy=TRUE.


gawdis(tall, w.type = "equal", groups = c(1, 2, 3, 3, 3), fuzzy = TRUE)
#> [1] "Running w.type=equal on groups=c(1,2,3,3,3)"
#> [1] "Running w.type=equal on groups=c(1)"
#> [1] "Running w.type=equal on groups=c(2)"
#> [1] "Running w.type=equal on groups=c(3,3,3)"
#>           sp1       sp2       sp3       sp4       sp5       sp6
#> sp2 0.2777778                                                  
#> sp3 0.5555556 0.6111111                                        
#> sp4 0.3888889 0.3333333 0.5000000                              
#> sp5 0.7333333 0.6555556 0.1777778 0.5000000                    
#> sp6        NA        NA        NA        NA        NA          
#> sp7 0.6666667 0.8333333 0.3333333 0.7222222 0.2888889        NA

Tutorial (ecología - funcional): análisis por rasgos (Parte 1)

Pablo Jorgensen, basado en Bello et al. (2021). Handbook of Trait-Based Ecology: from Theory to R Tools. Cambridge University Press

2023-05-29