Tutorial (ecología - diversidad): paquetes vegan e iNEXT

Introducción

En ecología, la diversidad de un ensamblaje es una medida que cuantifica el número de especies presentes en un área. En particular, el número de especies de un ensamblaje local es un índice intuitivo y natural de la estructura de la comunidad. En ecología aplicada y biología de la conservación, el número de especies que permanecen en una comunidad representa además un indicador valioso en la lucha por preservar y restaurar las comunidades perturbadas. Sin embargo, a pesar de nuestra familiaridad con la riqueza de especies, es una variable sorprendentemente difícil de medir. Casi sin excepción, la riqueza de especies no puede medirse con precisión ni estimarse directamente mediante la observación, porque el número de especies observado es un estimador sesgado hacia abajo de la riqueza de especies completa (total) de un ensamblaje (i.e. la subestima). Muchos estudios recientes siguen ignorando algunos de los problemas fundamentales de muestreo y medición que pueden comprometer la estimación precisa de la riqueza de especies (Gotelli y Colwell 2011).

En este tutorial veremos cómo solucionar los problemas asociados con la estimación del número de especies (\(S\)) para la caracterización de hábitats locales (diversidad \(\alpha\)), regiones (diversidad \(\gamma\)) y gradientes ambientales (diversidad \(\beta\)). Analizaremos modelos de muestreo para la riqueza de especies que tienen en cuenta el sesgo de submuestreo, ajustando o controlando las diferencias en el número de individuos y el número de muestras tomadas (rarefacción), así como los modelos que utilizan distribuciones de abundancia o incidencia (presencia/ausencia) para estimar el número de especies no detectadas (estimadores de riqueza asintótica). Los métodos de estimación de la riqueza de especies que veremos pueden aplicarse a conjuntos de organismos que han sido identificados por su genotipo, a especies o a algún rango taxonómico superior, como el género o la familia. Escribiremos “especies” para simplificar.

Dado que estamos hablando de la estimación de la riqueza de especies, asumimos que se han tomado una o más muestras (mediante recolección u observación) de uno o más hábitats o regiones para algún grupo o grupos de organismos o ensamblajes (por ejemplo, especies de invertebrados bentónicos, de peces, de macroalgas, etc). Distinguimos dos tipos de datos utilizados en los estudios de riqueza: (1) los datos de incidencia, en los que cada especie detectada en una muestra de un conjunto se anota simplemente como presente (o ausente), y (2) los datos de abundancia, en los que la abundancia de cada especie se cuenta dentro de cada muestra. Por supuesto, los datos de abundancia siempre pueden convertirse en datos de incidencia (pero no a la inversa).

Por su naturaleza, los datos de muestreo sólo documentan la presencia verificada de especies en las muestras. La ausencia de una especie concreta en una muestra puede representar una ausencia real (la especie no está presente en el conjunto) o una ausencia falsa (la especie está presente, pero no se detectó en la muestra). Los métodos de estimación de la riqueza para los datos de abundancia suponen que los organismos pueden ser muestreados e identificados como individuos distintos. En el caso de los organismos clonales y coloniales, como muchas especies de algas y corales, los individuos no siempre pueden separarse o contarse, pero los métodos diseñados para los datos de incidencia pueden utilizarse, no obstante, si la presencia de las especies se registra en cuadrículas o muestras estandarizadas.

La utilidad de los métodos de inferencia estadística que veremos dependen en gran medida de la suposición biológica de que la comunidad es “cerrada”, con un número total de especies invariable y una distribución de abundancia de especies constante. En una metacomunidad abierta (en la que el ensamblaje cambia de tamaño y composición a lo largo del tiempo), puede que no sea posible hacer inferencias válidas sobre la estructura de la comunidad a partir de una muestra instantánea en un punto del tiempo. Pocas comunidades reales, si es que hay alguna, son completamente “cerradas”. Pero muchas están lo suficientemente circunscritas como para que puedan utilizarse los estimadores de riqueza (con precaución y advertencias). Para todos los métodos y métricas que analizamos partimos del supuesto estadístico estrechamente relacionado de que el muestreo es con reemplazo. En términos de recopilación de datos de inventario de la naturaleza este supuesto significa que, o bien se registran los individuos sin eliminarse del conjunto (por ejemplo, el censo de corales en una parcela) o, si se eliminan, las proporciones restantes no cambian por el muestreo.

Consideremos un gráfico en el que el eje x es el número de individuos muestreados y el eje y es el número acumulado de especies registradas. Imaginemos que se toma una muestra de un individuo por vez en el campo, al azar. A medida que se muestrean más individuos, el número total de especies registradas en la muestra aumenta, y se genera una curva de acumulación de especies. Por supuesto, el primer individuo extraído representará exactamente una especie nueva en la muestra, por lo que todas las curvas de acumulación de especies basadas en organismos individuales se originan en el punto [1,1]. La probabilidad de extraer una nueva especie dependerá tanto del número total de especies del ensamblaje como de sus abundancias relativas. Cuantas más especies haya en el ensamblaje y cuanto más uniforme sea la distribución de la abundancia de especies, más rápidamente subirá esta curva. Por el contrario, si la distribución de la abundancia de especies es muy desigual (unas pocas especies comunes y muchas raras, por ejemplo), la curva subirá más lentamente incluso al principio. Esto ocurre porque la mayoría de los individuos muestreados representarán especies más comunes que ya se han añadido a la muestra, en lugar de otras más raras que aún no se han detectado. Independientemente de la distribución de la abundancia de las especies, esta curva aumenta de forma monótona con una pendiente desacelerada. Para una muestra determinada, diferentes realizaciones estocásticas del orden en que los individuos de la muestra se añaden al gráfico producirán curvas de acumulación de especies que difieren ligeramente entre sí.

Rarefacción

En principio trabajaremos con datos de conteo de árboles (que podrían ser de mangle) en parcelas de 1 hectárea en una isla de Panamá. Para la mayoría de los análisis usarmeos el paquete de R “vegan” (Oksanen et al. 2022), que es uno de los más usados en Ecología de Comunidades, y que además de proveer varias funciones para el análisis de la diversidad también ofrece una gran variadad de análisis de ordenamiento y otros métodos multivariados. Los datos del ejemplo están incluídos eb vegan, los leemos a continuación.

EJEMPLO:

Leemos los datos BCI y, utilizando el método de rarefacción, comparamos el número de especies \(S\) de BCI usando las funciones rarefy y rarecurve del paquete vegan (Fig. 1 der.).

rm(list=ls())  # limpia variables de la memoria.

# install.packages(vegan)
library(vegan)
data(BCI) # Carga los conjuntos de datos especificados.
# str(BCI)  # Estructura de datos.
# View(BCI)  # Si necesitan visualizar todo el dataframe.
?BCI      # 50 parcelas (filas) x 225 especies (columnas).
# Podemos ver cómo varía el número de individuos x especies y x parcela:
indXmuestra <- rowSums(BCI) # ?rowSums()
str(indXmuestra)
#>  Named num [1:50] 448 435 463 508 505 412 416 431 409 483 ...
#>  - attr(*, "names")= chr [1:50] "1" "2" "3" "4" ...
# head(indXmuestra)
summary(indXmuestra)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   340.0   409.0   428.0   429.1   443.5   601.0

Srar <- rarefy(BCI, min(indXmuestra))
# Srar

par(mfrow=c(1,2)) # Para visualizar 2 gráficos en una fila (2 columnas).
boxplot(indXmuestra, data=indXmuestra,
        col=c("blue"),
        main="Box-Plot",
        xlab="", ylab="Número de árboles") # frame.plot=TRUE,axes=FALSE

rarecurve(BCI, col = "blue")

Figura 1. Variabilidad del número de organismos por muestra (izq.), y curvas de rarefacción del número de especies de árboles en las parcelas boscosas de BCI (identificadas por los números enmarcados)

# Produce curvas de rarefacción. Los cuadrados marcan la 
# identidad de los sitios analizados. Para "enrarefaccionar" un sitio,
# siga la curva de rarefacción hasta que la curva se corresponda
# con el sitio de menor abundancia para obtener la riqueza
# de especies enrarecida.

Observamos que el 50% de las parcelas varían entre 409 y 443 árboles, y que la parcela de menor abundancia tuvo 340 árboles y la de mayor abundancia 601 árboles (figura 1, izq.). Entonces querremos comparar las parcelas por su riqueza estandarizando por el menor número de árboles observados en una parcela (esto es, 340 árboles). Esto es deseable ya que cuanto mayor sea el número de árboles contados en una parcela, más especies esperaré encontrar.

OBSERVACIÓN: La curva de acumulación de especies que describimos en la Introducción sugiere una forma intuitiva de comparar la riqueza de dos muestras que difieren en el número de individuos muestreados. Supongamos que una de las dos muestras tiene \(N\) individuos y \(S\) especies, y la otra tiene \(n\) individuos y \(s\) especies. Las muestras difieren en el número de individuos presentes (\(N>n\)) y suelen diferir en el número de especies presentes (normalmente \(S>s\)). En el procedimiento llamado “rarefacción”, extraemos al azar \(n^{\ast}\) individuos submuestreando sin reemplazo de la mayor de las dos muestras originales (donde \(n^{\ast} = n\), el tamaño de la muestra original más pequeña). Calculando el número medio de especies, \(\bar{s}^{\ast}\), entre submuestras repetidas de \(n^{\ast}\) individuos se estima \(E(s^\ast|n^{\ast})\), el número esperado de especies en una submuestra aleatoria de \(n^{\ast}\) individuos de la muestra original más grande. En otras palabras, la rarefacción comparará la riqueza de especies \(\hat{S}\) esperada para una muestra pequeña de \(n\) individuos extraída al azar del gran conjunto de \(N\) individuos. Utiliza el menor número de individuos observados en una muestra para interpolar el número esperado si todas las demás muestras sólo tuvieran ese número de individuos. Las curvas de rarefacción suelen crecer rápidamente al principio, a medida que se encuentran las especies más comunes, pero se estabilizan cuando sólo quedan por muestrear las especies más raras.

Se usó la función “rarefy” (vegan) para comparar, por el método de rarefacción, la riqueza específica esperada de las parcelas boscosas en sub-muestras aleatorias del tamaño especificado en el argumento “sample”, que es el segundo argumento especificado en el uso de la función. El argumento “sample” es igual al menor número de individuos observados en una parcela (Parcela 23, con 340 individuos), y sirve para estimar el número de especies esperado para ese número de individuos en las demás parcelas boscosas. La idea es interpolar la abundancia de las otras parcelas a la de menor abundancia para estimar \(\hat{S}\) con base a la misma cantidad de individuos. “Srar” es el objeto resultante de aplicar la función “rarefy” a los datos. “Srar” contiene el número esperado de especies para cada una de las 50 parcelas de BCI en muestras de 340 individuos. O sea, cuántas especies habría en cada parcela si todas tuvieran 340 individuos.

EJEMPLO:

Utilizando el método de rarefacción comparamos el número de especies \(S\) entre los sitios 19 y 35 de BCI usando la función “rarefy” del paquete vegan. Usamos “rarecurve” para graficar las curvas de rarefacción de ambos sitios. Interpretamos los resultados y concluímos cuál de los dos sitios es más diverso.

La función “rarefy” del paquete vegan se usa para comparar la \(\hat{S}\) estimada de las parcelas con un mismo número de individuos por interpolación. La abundancia de la parcela con menor número de individuos, entre las dos parcelas comparadas, es la 19 (433 árboles de 109 especies diferentes). En la parcela 35 se contabilizaron 601 árboles de 83 especies distintas. El número de especies estimado en cada parcela para una abundancis de 433 árboles fue: 109 para la parcela 19 (obviamente coincide con el número observado de especies), y 70 para la parcela 35 (suponiendo 433 árboles). Con base en las estimaciomnes de riqueza (\(\hat{S}\)) para un mismo número de árboles (433 individuos), la parcela 19 es más rica en especies que la 35 (109 vs. 70, respectivamente). Las diferencias de riqueza específica de ambas parcelas se ilustran mediante las curvas de rarefacción generadas con la función rarecurve (figura 2).

BCI.sel <- BCI[c("19", "35"), ]
indXBCI.sel <- rowSums(BCI.sel)
indXBCI.sel
#>  19  35 
#> 433 601

# Para contar el número de especies observado en cada parcela:
Sobs <- rowSums(BCI.sel != 0)
Sobs
#>  19  35 
#> 109  83

Srar <- rarefy(BCI.sel, min(indXBCI.sel))
Srar
#>        19        35 
#> 109.00000  69.71908 
#> attr(,"Subsample")
#> [1] 433

rarecurve(BCI.sel, col = "blue")

Figura 2. Curvas de rarefacción para las parcelas 19 y 35

OBSERVACIÓN: Lamentablemente las funciones del paquete vegan usadas no permiten la estimación de la variabilidad del número esperado de especies para \(n\) individuos. Sin embargo, con el paquete iNEXT es posible ajustar un intervalo de confianza paramétrico del 95\(\%\) (u otro) mediante bootstrap (siguiendo un procedimiento similar al descrito para los reordenamientos aleatorios de individuos). Así, puede verificarse si el intervalo de confianza del 95\(\%\) de la riqueza media de especies observada de la muestra completa más pequeña se solapa con el intervalo de confianza del 95\(\%\) de la riqueza de especies esperada basada en submuestras aleatorias de tamaño \(n\) de la muestra más grande. Si se solapan, entonces la hipótesis de que la riqueza de la muestra más pequeña (basada en todos los \(n\) individuos) no difiere de la riqueza de una submuestra de tamaño \(n^{\ast}\) de la muestra más grande no puede rechazarse (no rechazo la hipótesis nula a un \(p \leq 0,05\)). Si no se rechaza esta hipótesis nula y las muestras originales (no rarefaccionada) difieren en la densidad de especies, entonces esta diferencia en la densidad de especies está determinada por un número diferente de individuos entre las dos muestras. Alternativamente, si \(s\) no está contenido dentro del intervalo de confianza de \(s^{\ast}\), las dos muestras difieren en la riqueza de especies de manera que no puede explicarse enteramente por las diferencias en la abundancia y/o el esfuerzo de muestreo (a un \(p \leq 0,05\)).

EJEMPLO:

Utilizando la función iNEXT del paquete homónimo comparamos las curvas de rarefacción de (la transpuesta) BCI.sel especificando el argumento datatype = ‘abundance’. Con “plot()” reproducimos gráficamente las curvas de rarefacción y sus IC para ambos sitios. Nota: q0 se refiere al índice de riqueza específica \(S\) en la serie de índices de diversidad de Hill. Gotelli (1996) argumentó varias razones por las que es preferible usar \(S\), en lugar de los otros índices de la serie de Hill (Simpson, Shannon-Wiener), para describir a la diversidad de especies de un ensamble o comunidad.

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

BCI.sel <- BCI[c("19", "35"), ]
DataInfo(t(BCI.sel), datatype = 'abundance')

par(mfrow = c(1, 2))
# Para iNEXT debo transponer al data.frame:
D_abund.BCI.sel <- iNEXT(t(BCI.sel), datatype = 'abundance')
plot (D_abund.BCI.sel)

# Para estimar la diversidad S de los sitios respecto al menor número de individuos
# observados en un sitio (433 individuos en el sitio 19), podemos utilizar 
D_abund.BCI.min <- iNEXT(t(BCI.sel), datatype = 'abundance', q=0, endpoint = 433)
plot (D_abund.BCI.min)

Figura 3. Curvas de rarefacción para las parcelas 19 y 35. Se muestran los intervalos de confianza de las estimaciones del número de especies. Izq.: Las partes de las curvas segmentadas muestran las extrapolaciones hechas por la función. Der.: Las estimaciones se cortan en el número de individuos correspondiente a la parcela con menos árboles cuantificados (parcela 19)

Los números estimados se calculan, de forma más eficiente mediante la función estimateD, con casi los mismos argumentos que iNEXT:

# Los números estimados se calculan, de forma más eficiente, mediante
# la función estimateD:  
D_abund.BCI.est <- estimateD(t(BCI.sel), datatype = 'abundance', q=0)
D_abund.BCI.est

?estimateD

Interpreten los resultados de los dos gráficos de la figura 3 ¿cuál es la diferencia entre ambas?. La principal diferencia entre ambos gráficos realizados con la función plot sobre el objeto generado con iNEXT es que el de la izquierda muestra la totalidad de las curvas de ambas parcelas estudiadas. En esa figura (figura 3, izq.) las curvas segmentadas muestran las extrapolaciones hechas por la función para estimar el número de especies con variaciones del número de individuos más elevados que el número de individuos observados por parcela (mostrados con símbolos en la figura). En cambioo, en la figura 3 (derecha) las estimaciones del número de especies se cortan en el número de individuos correspondiente a la parcela con menos árboles cuantificados (parcela 19, con 433 árboles). Más allá de esas diferencias, ambos gráficos muestran diferencias evidentes en el número de especies estimado en ambas parcelas, a lo largo de casi la totalidad del esfuerzo de muestreo (definido por el número de inividuos). Notar que los intervalos de confianza de cada parcela no se solapan.

¿Qué podemos concuir teniendo en cuenta además la tabla de la salida con estimateD?. Concluyo que la diversidad de la parcela 19 es mayor que la diversidad de la 35. El número de especies estimado (\(\hat{S}\)) por extrapolación con estimateD (los “qD” en la tabla de salida de estimateD) indicó 145 especies para la parcela 19, contra 101 especies en la parcela 35, considerando una muestra de 866 indivíduos en ambas parcelas. Notar que 866 corresponde al número de individuos en el extremo derecho de la curva con menor número de individuos observados en una parcela (parcela 19, figura 3 izquierda). Las diferencias son significativas porque los intervalos de confianza de ambas estimaciones de \(\hat{S}\) no se solapan, como puede verse en la figura 3 (izq.) y, principalmente, en la salida de estimateD (valores de “qD.LCL” y “qD.UCL”, que corresponden a los intervalos de confianza inferior y superior, respectivamente).

¿Qué valor de tamaño de muestra utiliza como referencia para la comparación del número de especies?. La estimación del número de especies (\(\hat{S}\)) para las parcelas 19 y 35 está basada en una muestra de 866 indivíduos (árboles), por lo explicado en la pregunta anterior.

EJEMPLO:

Otro conjunto de datos que incluye el paquete vegan es “dune”, que cuantifica la cobertura de la vegetación de dunas en una región costera de Países Bajos. Para 20 sitios se proporcionan los valores de clase de cobertura de 30 especies de plantas, pero - para los fines de los ejercicios a continuación - supondremos que los números indicando valores de clase de cobertura por especie vegetal representan conteos de plantas por especies.

Leemos los datos “dune” y describimos la comunidad vegetal de dunas de los sitios estudiados en función de su riqueza específica. Las curvas de rarefacción del ensamblaje vegetal muestran una alta variabilidad en el esfuerzo de muestreo (por número de individuos) entre los sitios o muestras (sitios 1 a 20) (figura 4). Por ejemplo, mientras que en el sitio 6 se contabilizaron más de 40 individuos (plantas), en el sitio 17 se contabilizaron menos de la mitad de plantas (figura 4). Algunas curvas con abundancias intermedias de individuos (plantas) parecen cruzarse, limitando el análisis de rarefacción. Sin embargo varios sitios se diferencian claramente por su riqueza de especies, observada y estimada por las curvas, variando entre el de menor diversidad (sitio 1) y el sitio más diverso en especies vegetales (sitio 5) (figura 4). Para discernir estas diferencias es útil la salida de la función “rarefy” del paquete vegan, que permite comparar la \(\hat{S}\) estimada de las parcelas o sitios con base en un mismo número de individuos. Con base en 15 individuos por sitio, se comprueba que el sitio 1 es el de menor riqueza específica (\(\hat{S} \approx 5\) especies), aunque no el de menor abundancia, y que el sitio 5 es el más rico en especies (\(\hat{S} \approx 10\) especies).

data(dune)
?dune
str(dune)
#> 'data.frame':    20 obs. of  30 variables:
#>  $ Achimill: num  1 3 0 0 2 2 2 0 0 4 ...
#>  $ Agrostol: num  0 0 4 8 0 0 0 4 3 0 ...
#>  $ Airaprae: num  0 0 0 0 0 0 0 0 0 0 ...
#>  $ Alopgeni: num  0 2 7 2 0 0 0 5 3 0 ...
#>  $ Anthodor: num  0 0 0 0 4 3 2 0 0 4 ...
#>  $ Bellpere: num  0 3 2 2 2 0 0 0 0 2 ...
#>  $ Bromhord: num  0 4 0 3 2 0 2 0 0 4 ...
#>  $ Chenalbu: num  0 0 0 0 0 0 0 0 0 0 ...
#>  $ Cirsarve: num  0 0 0 2 0 0 0 0 0 0 ...
#>  $ Comapalu: num  0 0 0 0 0 0 0 0 0 0 ...
#>  $ Eleopalu: num  0 0 0 0 0 0 0 4 0 0 ...
#>  $ Elymrepe: num  4 4 4 4 4 0 0 0 6 0 ...
#>  $ Empenigr: num  0 0 0 0 0 0 0 0 0 0 ...
#>  $ Hyporadi: num  0 0 0 0 0 0 0 0 0 0 ...
#>  $ Juncarti: num  0 0 0 0 0 0 0 4 4 0 ...
#>  $ Juncbufo: num  0 0 0 0 0 0 2 0 4 0 ...
#>  $ Lolipere: num  7 5 6 5 2 6 6 4 2 6 ...
#>  $ Planlanc: num  0 0 0 0 5 5 5 0 0 3 ...
#>  $ Poaprat : num  4 4 5 4 2 3 4 4 4 4 ...
#>  $ Poatriv : num  2 7 6 5 6 4 5 4 5 4 ...
#>  $ Ranuflam: num  0 0 0 0 0 0 0 2 0 0 ...
#>  $ Rumeacet: num  0 0 0 0 5 6 3 0 2 0 ...
#>  $ Sagiproc: num  0 0 0 5 0 0 0 2 2 0 ...
#>  $ Salirepe: num  0 0 0 0 0 0 0 0 0 0 ...
#>  $ Scorautu: num  0 5 2 2 3 3 3 3 2 3 ...
#>  $ Trifprat: num  0 0 0 0 2 5 2 0 0 0 ...
#>  $ Trifrepe: num  0 5 2 1 2 5 2 2 3 6 ...
#>  $ Vicilath: num  0 0 0 0 0 0 0 0 0 1 ...
#>  $ Bracruta: num  0 0 2 2 2 6 2 2 2 2 ...
#>  $ Callcusp: num  0 0 0 0 0 0 0 0 0 0 ...
dune.indXmue <- rowSums(dune)  # ?rowSums()
dune.indXmue
#>  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 
#> 18 42 40 45 43 48 40 40 42 43 32 35 33 24 23 33 15 27 31 31
summary(dune.indXmue)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   15.00   30.00   34.00   34.25   42.00   48.00

Srar <- rarefy(dune, min(dune.indXmue))
Srar
#>        1        2        3        4        5        6        7        8 
#> 4.813725 8.294747 7.985897 9.105309 9.787891 8.688685 9.496061 9.317495 
#>        9       10       11       12       13       14       15       16 
#> 9.575586 9.027570 7.734350 7.719151 7.782660 6.572498 7.271063 6.961173 
#>       17       18       19       20 
#> 7.000000 7.752880 7.891114 7.334569 
#> attr(,"Subsample")
#> [1] 15

rarecurve(dune, col = "blue")

Figura 4. Curvas de rarefacción para el número de especies vegetales en dunas (los números enmarcados identifican a la muestra

Utilizando las funciones y métodos apropiados comparamoa formalmente la riqueza de especies de los sitios 4, 12 y 17. ¿Varía la diversidad de los sitios?. La mejor opción para comparar formalmente (inferencialmente, a través de pruebas de hipótesis y la estimación de parámetros) la diferencia de diversidad entre los tres sitios la ofrece el paquete iNEXT. La función iNEXT genera un objeto que permite visualizar las curvas de rarefacción y la extrapolación del número de especies estimado para abundancias superiores al número de individuos observados en cada sitio o muestra. Esto representa una ventaja importante sobre la función rarefy de vegan. Vegan no ofrece la posibilidad de comparar la riqueza de especies por métodos de extrapolación con mediddas de error, sólo permite la interpolación (permite la extrapolación, pero con error estándar = 0). Los gráficos generados a partir de la lectura del objeto de iNEXT ilustran los intervalos de confianza alrededor de las estimaciones de \(\hat{S}\), lo que permite una evaluación rápida de las diferencias de diversidad \(\alpha\) entre sitios. Sin embargo, los datos para la comparación formal son generados a partir de la aplicación de la función ëstimateD del paquete iNEXT sobre la matriz, o “data frame”, de los datos comunitarios. La aplicación de esta función genera una tabla con información de la estimación de \(\hat{S}\), su intervalo de confianza alrededor de esa estimación, y el tamaño muestral sobre el que se realizaron las estimaciones (o sea, el número de individuos de referencia para las estimaciones de riqueza específica), entre otras estimaciones. Este tamaño muestral, como mencioné antes, puede ser una extrapolación del número de individuos observado para las parcelas o sitios analizados. El método de estimación de los parámetros deseados es también reportado (extrapolación o rarefacción).

Aplicando la función iNEXT sobre la matriz de datos de la vegetación de dunas de los sitios, 4, 12 y 17, y graficando el objeto devuelto, generé la figura 5 (izq.); pero también obtuve información básica del total de sitios con la función “DataInfo”. Repasando, puedo verificar que en el sitio 4 se observaron 13 especies (de 45 plantas en la muestra), en el sitio 12 se observaron 9 especies (de 35 plantas en las muestra) y en el sitio 17 se contaron 7 especies (en 15 plantas individuales que conformaron la muestra). Más allá de la variación en el número de especies y el número de plantas entre muestras, puede verse que el sitio 4 parece ser significativamente más rico en especies que los otros dos sitios, mientras que los sitios 12 y 17 no se diferenciarían en su riqueza específica (figura 5). Nótese que para la evaluación de las diferencias en la roqueza de especies entre los sitios hice el corte imaginario en el extremo derecho de la curva azul (sitio 17) - la de menor número de individuos por muestra entre los sitios comparados - correspondiente a un tamaño de nuestra de alrededor de 30 individuos (figura 5, izquierda). No hice el corte en el número observado de especies para ese sitio (indicado por el símbolo azul), como es habitual con el método de rarefacción (función rarefy del paquete vegan) (figura 5, derecha).

Para comparar formalmente la \(\hat{S}\) en los tres sitios especificados observo las estimaciones con un mismo número de individuos, de la salida de estimateD. En la tabla de la salida de esa función del paquete iNEXT se confirman los resultados descritos a partir del análisis visual (figura 5), ya que el número de especies estimado para una muestra de 30 individuos (m) es: 12.1 especies (sitio 4), 9.0 (sitio 12) y 7.1 especies (sitio 17). Para ver si esas diferencias son significativas observo si se solopan los intervalos de confianza de las estimaciones anteriores de \(\hat{S}\) (qD de la tabla debajo). Los límites de los intervalos son definidos por qD.LCL y qD.UCL (límites inferior y superior, respectivamente). Se confirma también el patrón descrito a partir de la figura 5, el sitio 4 tiene más especies que cualquiera de los otros dos sitios (12 y 17), pero los sitios 12 y 17 no se diferencia entre sí por su riqueza específica.

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

data(dune)
DataInfo(t(dune), datatype = 'abundance')

dune.sel <- dune[c("4", "12", "17"), ]
indXdune.sel <- rowSums(dune.sel)
indXdune.sel
#>  4 12 17 
#> 45 35 15

par(mfrow = c(1, 2))

# Para iNEXT debo transponer al data.frame:
D_abund.dune.sel <- iNEXT(t(dune.sel), datatype = 'abundance')
plot(D_abund.dune.sel)

# Para estimar la diversidad S de los sitios respecto al menor número de individuos
# observados en un sitio (433 individuos en el sitio 19), podemos utilizar 
D_abund.dune.sel.min <- iNEXT(t(dune.sel), datatype = 'abundance', q=0,
                              endpoint = min(indXdune.sel))
plot (D_abund.dune.sel.min)

Figura 5. Curvas de rarefacción para las muestras 4, 12 y 17. Se muestran los intervalos de confianza de las estimaciones del número de especies. Izq.: Las partes de las curvas segmentadas muestran las extrapolaciones hechas por la función. Der.: Las estimaciones se cortan en el número de individuos correspondiente a la muestra con menos plantas cuantificados (muestra o cuadrante 17)

# Los números estimados se calculan, de forma más eficiente, mediante
# la función estimateD:  
D_abund.dune.sel.est <- estimateD(t(dune.sel), datatype = 'abundance', q=0)
D_abund.dune.sel.est

OBSERVACIÓN: Las curvas de rarefacción suelen considerarse una solución objetiva para comparar la riqueza de especies con diferentes tamaños de muestra (o individuos). Sin embargo, el ránking entre los sitios por su riqueza específica puede diferir al establecer el corte a diferentes tamaños de muestra de rarefacción, ya que las curvas de rarefacción pueden cruzarse.

Modelos de acumulación de especies

En el paquete vegan los modelos de acumulación de especies (SAC) están enfocados a la estimación de la riqueza \(\gamma\) (riqueza de grupos de sitios). Los modelos SAC son similares a la rarefacción, estudian la acumulación de especies cuando el número de sitios aumenta. Existen varios métodos alternativos, como la acumulación de sitios en el orden en que se encuentran, y la acumulación repetida en orden aleatorio. Hay tres modelos analíticos, pero el recomendado es el método exacto de Kindt. La función specaccum (vegan) encuentra las curvas de acumulación de especies para un determinado número de sitios o individuos muestreados (figura 6).

EJEMPLO:

sac <- specaccum(BCI)
plot(sac, ci.type="polygon", ci.col="yellow")

Figura 6. Acumulación de especies para la comunidad estudiada de árboles en las parcelas (sitios) de BCI

EJEMPLO:

Obtenemos la diversidad \(\gamma\) de la comunidad vegetal de las dunas vistas en el ejemplo anterior. La diversidad \(\gamma\) sirve para cuantificar el número de especies de toda una región. La función “specaccum”, del paquete vegan, permite obtener curvas de acumulación del número de especies para un cierto número de sitios y es ideal para la estimación de la riqueza \(\gamma\). La curva de la figura 7 muestra la acumulación de especies con el aumento del esfuerzo de muestreo (número de sitios). La curva es asintótica con \(\gamma = 30\) especies en 20 sitios.

data(dune)
sac.dune <- specaccum(dune)
plot(sac.dune, ci.type="polygon", ci.col="orange")

Figura 7. Acumulación de especies vegetales en las muestras de la comunidad de dunas

?specaccum

Un argumento importante de la función es “method”, que sirve para especificar los algoritmos usados en la estimación de la curva y sus medidas de error. En el ajercicio anterior usé el método por defecto, que es method=“random”, que obtiene la media de las SAC y su desviación estándar a partir de permutaciones aleatorias de los datos, o submuestreo sin reemplazo. Resulta de interés la posibilidad de “pesar” la influencia de ciertos sitios con el argumento “w”. Los pesos especificados en “w” no influyen en el orden de acumulación de los sitios, sino que sólo diferencian a los sitios en función del esfuerzo de muestreo. El esfuerzo de muestreo real se indica como elemento “Effort” o “Individuals” en la salida de la función. Para el método “random” ponderado, el esfuerzo se refiere al esfuerzo medio por sitio, o a la suma de los pesos por número de sitios.

Diversidad Beta (\(\beta\))

Los índices de diversidad básicos son índices de diversidad alfa (\(\alpha\)), la diversidad en un punto. La diversidad beta (\(\beta\)) debe estudiarse con respecto a gradientes, pero casi todo el mundo la entiende como una medida de heterogeneidad general: ¿cuántas especies más hay en una colección de sitios en comparación con un sitio promedio?. El índice de diversidad \(\beta\) más conocido se basa en la relación entre el número total de especies en una colección de sitios, o gama (\(\gamma\)) y la riqueza media por un sitio (\(\bar{\alpha}\):

\[\beta=\gamma/\bar{\alpha}-1\]

La sustracción de uno significa que \(\beta = 0\) cuando no hay exceso de especies o no hay heterogeneidad entre sitios. Para este índice, no se necesitan funciones específicas, pero se puede calcular fácilmente con la ayuda de la función de vegan specnumber:

EJEMPLO:

ncol(BCI)/mean(specnumber(BCI)) - 1
#> [1] 1.478519

El índice \(\beta\) así calculado es problemático porque \(\gamma\) aumenta con el número de sitios incluso cuando los sitios son todos subconjuntos de la misma comunidad. Una alternativa se construye haciendo una comparación por pares de sitios con el índice de disimilitud de Sorensen, que puede calcularse para todos los sitios utilizando la función “vegdist” (con datos binarios):

beta <- vegdist(BCI, binary=TRUE)
mean(beta)
#> [1] 0.3399075

Existen muchas otras definiciones para la diversidad beta. Todos los índices comúnmente utilizados pueden hallarse utilizando betadiver, que pueden ser referidos por el nombre del subíndice, o el número del índice:

betadiver(help=TRUE)
#> 1 "w" = (b+c)/(2*a+b+c)
#> 2 "-1" = (b+c)/(2*a+b+c)
#> 3 "c" = (b+c)/2
#> 4 "wb" = b+c
#> 5 "r" = 2*b*c/((a+b+c)^2-2*b*c)
#> 6 "I" = log(2*a+b+c) - 2*a*log(2)/(2*a+b+c) - ((a+b)*log(a+b) +
#> (a+c)*log(a+c)) / (2*a+b+c)
#> 7 "e" = exp(log(2*a+b+c) - 2*a*log(2)/(2*a+b+c) - ((a+b)*log(a+b) +
#> (a+c)*log(a+c)) / (2*a+b+c))-1
#> 8 "t" = (b+c)/(2*a+b+c)
#> 9 "me" = (b+c)/(2*a+b+c)
#> 10 "j" = a/(a+b+c)
#> 11 "sor" = 2*a/(2*a+b+c)
#> 12 "m" = (2*a+b+c)*(b+c)/(a+b+c)
#> 13 "-2" = pmin(b,c)/(pmax(b,c)+a)
#> 14 "co" = (a*c+a*b+2*b*c)/(2*(a+b)*(a+c))
#> 15 "cc" = (b+c)/(a+b+c)
#> 16 "g" = (b+c)/(a+b+c)
#> 17 "-3" = pmin(b,c)/(a+b+c)
#> 18 "l" = (b+c)/2
#> 19 "19" = 2*(b*c+1)/(a+b+c)/(a+b+c-1)
#> 20 "hk" = (b+c)/(2*a+b+c)
#> 21 "rlb" = a/(a+c)
#> 22 "sim" = pmin(b,c)/(pmin(b,c)+a)
#> 23 "gl" = 2*abs(b-c)/(2*a+b+c)
#> 24 "z" = (log(2)-log(2*a+b+c)+log(a+b+c))/log(2)
# ?betadiver()

Algunos de estos índices son duplicados, y muchos de ellos son índices de disimilitud bien conocidos. Uno de los índices más interesantes se basa en el modelo de especies\(\sim\)area de Arrhenius (\(SAR\)):

\[\hat{S} = c \cdot A^{z}\]

donde \(A\) es el tamaño del área del sitio o parche de hábitat; \(c\) es una constante derivada empíricamente que cuantifica el número de especies en una muestra por unidad de superficie dentro del sistema estudiado; y el exponente \(z\) es la pendiente de la relación especie-área en un gráfico log-log y, esto es lo importante, es una medida de la diversidad \(\beta\).
En las islas, normalmente \(z \approx 0,3\). Este tipo de islas puede considerarse como subconjuntos de una misma comunidad, lo que indica que realmente debemos hablar de diferencias de gradiente si \(z \gtrapprox 0,3\). Podemos encontrar el valor de z para un par de sitios o parcelas utilizando la función betadiver:

z <- betadiver(BCI, "z")
summary(z)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.2733  0.3895  0.4192  0.4213  0.4537  0.5906

El tamaño \(A\) y el parámetro \(c\) se cancelan, y el índice da la estimación \(z\) para cualquier par de sitios.

Obtenemos \(\beta\) para la comunidad vegetal de las dunas a partir de la ecuación de Arrhenius (\(SAR\)) ¿qué significa el \(\beta\) obtenido?. El índice \(\beta\) representa el recambio de especies a lo largo de un gradiente, es decir, nos expresa con qué tasa cambian las especies en relación con el cambio ambiental. Para la vegetación de dunas \(\beta = 0.62\), lo que indica que existe un recambio de especies fuerte en la región, probablemente en respuesta a algún gradiente ambiental marcado en el área del muestreo.

z <- betadiver(dune, "z")
summary(z)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.1155  0.4873  0.6374  0.6212  0.7655  1.0000

Banco de especies: Número de especies no vistas

Tamaño del banco para múltiples sitios

Los modelos de acumulación de especies indican que no todas las especies fueron observadas en algún sitio. Estas especies no contempladas también pertenecen al conjunto (o “pool”) de especies. Las funciones “specpool” y “estimateR” implementan algunos métodos para estimar el número de especies no observadas. La función “specpool” examina una colección de sitios, y “estimateR” trabaja con recuentos de individuos, y puede utilizarse con un solo sitio. Ambas funciones asumen que el número de especies no observadas está relacionado con el número de especies raras (especies avistadas sólo una o dos veces).

Las funciones basadas en la incidencia agrupan las especies por su número de apariciones. La proporción de muestreo \(i\) es una estimación de la frecuencia de la especie en la comunidad. Cuando las especies están presentes en la comunidad pero no en la muestra, \(i = 0\) es una subestimación obvia. Consecuentemente, para los valores \(i > 0\) la comunalidad de las especies está sobreestimada. Los modelos para el tamaño del pool estiman el número de especies que faltan en la muestra (\(f_0\)). La función specpool implementa modelos para estimar el número de especies que faltan (\(f_0\)). El método utiliza estimadores jackknife de primer y segundo orden y bootstrap. La idea en jackknife es que nos perdimos tantas especies como las que vimos una sola vez. La idea en bootstrap es que si repetimos el muestreo, (con reemplazo) de los mismos datos, nos perdemos tantas especies como las que nos perdimos originalmente. Los estimadores de la varianza sólo se refieren al número estimado de especies que faltan \(\hat{f_0}\). La varianza del jackknife de primer orden se basa en el número de “singles” (especies que aparecen una sola vez en los datos) en los gráficos de la muestra.

Los valores de riqueza extrapolados para todos los datos de BCI son:

EJEMPLO:

specpool(BCI)

Si la estimación del tamaño del conjunto funciona realmente, deberíamos obtener los mismos valores de riqueza estimada si tomamos un subconjunto aleatorio de la mitad de las parcelas (pero esto no suele ser cierto):

s <- sample(nrow(BCI), 25)
s
#>  [1] 15 22 35 40 18 41 33 24 10 31 43 14  1 28 19 47 39 46 20  2 13 23 16 37  5
specpool(BCI[s,])

¿Cuántas especies habría en realidad en la comunidad de dunas?. Expliquen su respuesta. La función “specpool” reporta cuatro estimadores del número real de especies que habría en la comunidad de dunas (teniendo en cuenta las especies “no vistas” o muestreadas). Los estimadores dan entre 32 y 33 especies aproximadamente. O sea, 2 a 3 especies más de las observadas en la comunidad. Esto significa que hay 2 ó 3 especies que no fueron muestreadas, pero que estarían en el pool de especies regional.

specpool(dune)

Tamaño del banco para un sitio

La función specpool necesita una colección de sitios, pero hay algunos métodos que estiman el número de especies no observadas para cada sitio individualmente. Estas funciones necesitan recuentos de individuos, y las especies observadas sólo una o dos veces (u otras especies raras) ocupan el lugar de las especies con frecuencias bajas. La función estimateR implementa dos de estos métodos:

EJEMPLO:

k=25
estimateR(BCI[k,])
#>                  25
#> S.obs    105.000000
#> S.chao1  142.142857
#> se.chao1  15.529860
#> S.ACE    148.773956
#> se.ACE     6.536657

El método de Chao es similar al modelo con corrección de sesgo. El estimador de la varianza utiliza una solución iterativa. El tamaño del “pool” se estima por separado para cada sitio (si la entrada es un data frame, se analizará cada sitio). Si el modelo de abundancia log-normal es apropiado para los datos, se puede utilizar para estimar el tamaño del pool. El modelo log-normal tiene un número definido de especies que se puede encontrar integrando la log-normal. La función veiledspec estima esta integral a partir de un modelo ajustado con prestondistr (o prestonfit). Utiliza prestonfit si se dan los datos brutos del sitio. El modelo log-normal puede ser deficiente, pero asuminos que describe bien nuestros datos de abundancia:

veiledspec(prestondistr(BCI[k,]))
#> Extrapolated     Observed       Veiled 
#>    118.48959    105.00000     13.48959
veiledspec(BCI[k,])
#> Extrapolated     Observed       Veiled 
#>    125.04651    105.00000     20.04651

¿Cuántas especies habría en uno de los sitios de la comunidad de dunas?. La función “estimateR” indica que para la muestra - o sitio - 10 habría entre 0 y 2 especies más, no contabilizadas en el sitio por el muestreo.

# Con estimateR:  
k=10
estimateR(dune[k,])
#>                  10
#> S.obs    12.0000000
#> S.chao1  12.0000000
#> se.chao1  0.1595712
#> S.ACE    12.2857143
#> se.ACE    0.6341603

# Con veiledspec:  
veiledspec(prestondistr(dune[k,]))
#> Extrapolated     Observed       Veiled 
#> 1.200075e+01 1.200000e+01 7.468542e-04
veiledspec(dune[k,])
#> Extrapolated     Observed       Veiled 
#>    13.973705    12.000000     1.973705

Probabilidad de pertenencia al pool

El suavizado de Beals se sugirió originalmente como una herramienta de regularización de datos para la ordenación. Regulariza demasiado los datos, pero se ha sugerido como método para estimar cuáles de las especies que faltan podrían darse en un sitio, o qué sitios son adecuados para una especie. La probabilidad de cada especie en cada sitio se evalúa a partir de otras especies que se dan en el sitio. La función beals implementa el suavizado de Beals (McCune, 1987; De Caceres y Legendre, 2008).

EJEMPLO:

smo <- beals(BCI)
# smo

Podemos ver cómo se relacionan la probabilidad estimada de aparición y el número observado de árboles en una de las especies más conocidas (figura 8). Estudiamos sólo una especie, y para evitar un razonamiento circular no incluimos la especie objetivo en el suavizado (Fig. 8).

j <- which(colnames(BCI) == "Ceiba.pentandra")
plot(beals(BCI, species=j, include=FALSE), BCI[,j], ylab="Occurrence",
     main="Ceiba pentandra", xlab="Probability of occurrence") 

Figura 8. Relación entre la probabilidad estimada de aparición y el número observado de árboles para una de las especies más conocidas