iNEXT: un paquete para estimar biodiversidad

La cuantificación o medición de la biodiversidad es un objetivo muy recurrente en estudios de ecología, monitoreo biológico y proyectos de gestión ambiental. Sin embargo, cuantificar la biodiversidad de un sitio puede llegar a ser muy difícil o casi imposible de lograr. Por lo tanto, la medida más simple y más utilizada para cuantificarla es la riqueza de especies (número de especies presentes en un conjunto) (Chao et al, 2014).

La riqueza de especies puede ser un índice problemático debido a la intensidad del muestreo y a la distribucion de la abundancia de las especies. Asimismo, es un índice muy sensitivo al tamaño de la muestra, ya que la mayoria de las especies pertenecientes a un ensamble son “raras” o poco detectables, provocando que la muestra sea incompleta. Por lo tanto, desde el punto de vista estadístico, la riqueza de especies es muy dificil de estimar con precisión a partir de una muestra finita (Carmona-Galindo y Carmona, 2013).

iNEXT es un paquete que permite estimar la diversidad de especies utilizando Hill numbers de orden q: riqueza de especies (q=0), diversidad de Shannon (q=1, la exponencial de la entropía de Shannon) y diversidad de Simpson (q=2, la inversa de la concentración de Simpson) (Espinosa, 2003). Para cada medida de diversidad, iNEXT utiliza la muestra observada de datos de abundancia o incidencia (presencia/ausencia) para calcular estimaciones de diversidad y los intervalos de confianza asociados del 95% (predeterminados), así como trazar 2 tipos de curvas de rarefacción y extrapolación ( R/E):

  1. Basadas en el tamaño de la muestra: número de individuos en una #muestra / número de unidades de muestreo para datos de incidencia.

  2. Basadas en cobertura de muestreo: traza las estimaciones de #diversidad con respecto a la cobertura de la muestra.

install.packages(“iNEXT”) ## Paquetes necesarios:

library(devtools)
## Warning: package 'devtools' was built under R version 4.3.2
## Loading required package: usethis
## Warning: package 'usethis' was built under R version 4.3.2
library(iNEXT)
## Warning: package 'iNEXT' was built under R version 4.3.2
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.2
library(usethis)

El paquete de iNEXT cuenta con 4 bases de datos:
spider (lista de dos vectores) → se corre con datatype =“abundance” bird (formato data.frame) → se corre con datatype =“abundance” ant (lista de 5 vectores) → se corre con datatype =“incidence_freq” ciliates (lista de 3 matrices) → se corre con datatype =“incidence_raw”

Ejemplo cuando tenemos datos de incidencia:

La base de datos “ant” da la frecuencia de incidencia de especies de hormigas en cinco elevaciones/conjuntos (Longino y Colwell, 2011). La unidad de muestreo es una parcela de suelo de bosque de 1m x 1m, para las 5 elevaciones de: 50 m, 500 m, 1070 m, 1500 m y 2000 m. El número de especies observadas en cada elevación fue de 227, 241, 122, 56 y 14 respectivamente.

data(ant)
str(ant)
## List of 5
##  $ h50m  : num [1:228] 599 330 263 236 222 195 186 183 182 129 ...
##  $ h500m : num [1:242] 230 133 131 123 78 73 65 60 60 56 ...
##  $ h1070m: num [1:123] 150 99 96 80 74 68 60 54 46 45 ...
##  $ h1500m: num [1:57] 200 144 113 79 76 74 73 53 50 43 ...
##  $ h2000m: num [1:15] 200 80 59 34 23 19 15 13 8 8 ...

Ejemplo cuando tenemos datos de abundancia:

La base de datos “spider” consiste en la abundancia de especies de arañas colectadas en dos tratamientos (“Girdled” y “Logged”).

data(spider) 
str(spider) 
## List of 2
##  $ Girdled: num [1:26] 46 22 17 15 15 9 8 6 6 4 ...
##  $ Logged : num [1:37] 88 22 16 15 13 10 8 8 7 7 ...

¿Cuál sería nuestro número de especies esperadas si seguimos muestreando?

Este comando nos dice el número de especies obtenidas y esperadas si se sigue muestreando. Para base de datos de incidencia (tipo incidence_freq):

ChaoRichness(ant, datatype = "incidence_freq", conf = 0.95)
##        Observed Estimator Est_s.e. 95% Lower 95% Upper
## h50m        227   279.109   19.794   252.374   334.009
## h500m       241   314.810   23.259   281.384   375.902
## h1070m      122   146.337   12.082   131.699   183.063
## h1500m       56    77.019   16.353    61.458   136.949
## h2000m       14    14.249    0.726    14.013    18.716

Para base de datos de abundancia:

ChaoRichness(spider, datatype = "abundance", conf = 0.95)
##         Observed Estimator Est_s.e. 95% Lower 95% Upper
## Girdled       26    43.893   14.306    30.511    96.971
## Logged        37    61.403   18.532    43.502   128.583

¿Cómo obtener el índice de diversidad de Shannon?

Para base de datos de incidencia (tipo incidence_freq):

ChaoShannon(ant, datatype = "incidence_freq", transform = FALSE, conf = 0.95, B = 200)
##        Observed Estimator Est_s.e 95% Lower 95% Upper
## h50m      4.376     4.400   0.015     4.376     4.430
## h500m     4.622     4.678   0.020     4.638     4.718
## h1070m    4.085     4.128   0.026     4.085     4.178
## h1500m    3.283     3.315   0.026     3.283     3.366
## h2000m    2.067     2.092   0.051     2.067     2.193

Para base de datos de abundancia:

ChaoShannon(spider, datatype = "abundance", transform = FALSE, conf = 0.95, B = 200)
##         Observed Estimator Est_s.e 95% Lower 95% Upper
## Girdled    2.490     2.627   0.110     2.490     2.841
## Logged     2.669     2.793   0.107     2.669     3.004

Tengamos presente que los valores de diversidad de este índice de Shannon van de 0.5 a 5, los valores más comunes son entre 2 y 3, siendo los menores de 2 poco diversos y mayores a 3 diversos.

¿Cómo obtener el índice de diversidad de Simpson?

Para base de datos de incidencia (tipo incidence_freq):

ChaoSimpson(ant, datatype = "incidence_freq", transform = FALSE, conf = 0.95, B = 200)
##        Observed Estimator Est_s.e. 95% Lower 95% Upper
## h50m      0.980     0.980    0.000     0.980     0.981
## h500m     0.984     0.985    0.000     0.984     0.985
## h1070m    0.976     0.976    0.001     0.976     0.978
## h1500m    0.947     0.947    0.002     0.947     0.951
## h2000m    0.830     0.832    0.012     0.830     0.856

Para datos de abundancia:

ChaoSimpson(spider, datatype = "abundance", transform = FALSE, conf = 0.95, B = 200)
##         Observed Estimator Est_s.e. 95% Lower 95% Upper
## Girdled    0.872     0.878    0.017     0.872     0.911
## Logged     0.852     0.855    0.019     0.852     0.893

Mientras el valor de diversidad de simpson sea más cercano a 1 significa que hay dominancia de una especie sobre las demás, si el número es más cercano a 0 significa que no hay dominancia y es más diverso el sitio.

¿Cómo resumimos datos generales de iNEXT?

Para tener un resumen de los datos en un datatype = “abundance”:

DataInfo(spider, datatype="abundance")
##   Assemblage   n S.obs     SC f1 f2 f3 f4 f5 f6 f7 f8 f9 f10
## 1    Girdled 168    26 0.9289 12  4  0  1  0  2  0  1  1   0
## 2     Logged 252    37 0.9446 14  4  4  3  1  0  3  2  0   1

Donde: “n” es el número de individuos observados o colectados; “s.obs” es el número de especies encontradas; estimación de cobertura de muestra para la muestra de referencia (SC); los primeros diez recuentos de frecuencia (f1‐f10).

Para tener un resumen de los datos en un datatype= “incidence_freq” or “incidence_raw”:

DataInfo(ant, datatype="incidence_freq")
##   Assemblage   T    U S.obs     SC Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10
## 1       h50m 599 5976   227 0.9918 49 23 18 14  9 10  4  8  6   2
## 2      h500m 230 2943   241 0.9760 71 34 12 14  9 11  8  4  7   5
## 3     h1070m 150 1730   122 0.9839 28 16 13  3  1  3  6  1  1   1
## 4     h1500m 200 1170    56 0.9889 13  4  2  2  4  2  0  0  4   0
## 5     h2000m 200  271    14 0.9964  1  2  1  1  0  0  0  2  0   0

Donde: “T” es el tamaño de la muestra, “S.obs” la riqueza de especies observada, “U” es número total de incidencias, “SC” es estimación de cobertura de la muestra para la muestra de referencia y Q1‐Q10 son los primeros diez recuentos de frecuencia de incidencia.

Al aplicar el comando iNext este nos da 3 listas de salida:

iNEXT(spider, q=0, datatype="incidence_freq", conf = 0.95)
## Compare 2 assemblages with Hill number order q = 0.
## $class: iNEXT
## 
## $DataInfo: basic data information
##   Assemblage  T   U S.obs     SC Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10
## 1    Girdled 46 122    25 0.9031 12  4  0  1  0  2  0  1  1   0
## 2     Logged 88 164    36 0.9152 14  4  4  3  1  0  3  2  0   1
## 
## $iNextEst: diversity estimates with rarefied and extrapolated samples.
## $size_based (LCL and UCL are obtained for fixed size.)
## 
##    Assemblage   t        Method Order.q        qD    qD.LCL    qD.UCL        SC
## 1     Girdled   1   Rarefaction       0  2.652174  2.179985  3.124362 0.2480874
## 10    Girdled  23   Rarefaction       0 17.943782 14.667720 21.219843 0.8635463
## 20    Girdled  46      Observed       0 25.000000 19.861351 30.138649 0.9030753
## 30    Girdled  68 Extrapolation       0 29.867282 22.845586 36.888978 0.9298666
## 40    Girdled  92 Extrapolation       0 33.656399 24.148621 43.164178 0.9507233
## 41     Logged   1   Rarefaction       0  1.863636  1.574869  2.152404 0.1040090
## 50     Logged  44   Rarefaction       0 27.302142 23.468282 31.136003 0.8604931
## 60     Logged  88      Observed       0 36.000000 29.471735 42.528265 0.9151912
## 70     Logged 130 Extrapolation       0 41.822811 32.651149 50.994473 0.9355790
## 80     Logged 176 Extrapolation       0 46.607061 34.263739 58.950382 0.9523305
##        SC.LCL    SC.UCL
## 1  0.18968168 0.3064932
## 10 0.81635178 0.9107409
## 20 0.85789635 0.9482542
## 30 0.87793051 0.9818026
## 40 0.89829684 1.0000000
## 41 0.07548903 0.1325289
## 50 0.81669415 0.9042920
## 60 0.87308590 0.9572965
## 70 0.88956614 0.9815919
## 80 0.90580690 0.9988540
## 
## NOTE: The above output only shows five estimates for each assemblage; call iNEXT.object$iNextEst$size_based to view complete output.
## 
## $coverage_based (LCL and UCL are obtained for fixed coverage; interval length is wider due to varying size in bootstraps.)
## 
##    Assemblage        SC   t        Method Order.q        qD    qD.LCL    qD.UCL
## 1     Girdled 0.2480914   1   Rarefaction       0  2.652214  2.173472  3.130956
## 10    Girdled 0.8635464  23   Rarefaction       0 17.943787 10.202370 25.685203
## 20    Girdled 0.9030753  46      Observed       0 25.000000  8.972236 41.027764
## 30    Girdled 0.9298666  68 Extrapolation       0 29.867282  6.665179 53.069385
## 40    Girdled 0.9507233  92 Extrapolation       0 33.656399  4.070811 63.241987
## 41     Logged 0.1040095   1   Rarefaction       0  1.863645  1.594542  2.132748
## 50     Logged 0.8604930  44   Rarefaction       0 27.302138 18.134316 36.469961
## 60     Logged 0.9151912  88      Observed       0 36.000000 13.254542 58.745458
## 70     Logged 0.9355790 130 Extrapolation       0 41.822811 11.911993 71.733629
## 80     Logged 0.9523305 176 Extrapolation       0 46.607061 10.302620 82.911502
## 
## NOTE: The above output only shows five estimates for each assemblage; call iNEXT.object$iNextEst$coverage_based to view complete output.
## 
## $AsyEst: asymptotic diversity estimates along with related statistics.
##   Assemblage         Diversity Observed Estimator      s.e.       LCL      UCL
## 1    Girdled  Species richness 25.00000  42.60870 16.589647 25.000000 75.12381
## 2    Girdled Shannon diversity 13.74017  16.39728  1.751652 12.964102 19.83045
## 3    Girdled Simpson diversity 10.02965  10.69048  1.010397  8.710137 12.67082
## 4     Logged  Species richness 36.00000  60.22159 14.366841 36.000000 88.38008
## 5     Logged Shannon diversity 22.34412  26.90553  2.652376 21.706972 32.10410
## 6     Logged Simpson diversity 16.32039  17.91803  1.830347 14.330619 21.50545

$DataInfo $iNextEst para mostrar estimaciones de diversidad basadas en tamaño y cobertura junto con estadísticas relacionadas para una serie de muestras enrarecidas y extrapoladas (datos que se utilizan para la graficación). $AsyEst para mostrar estimaciones de diversidad asintótica junto con estadísticas relacionadas.

  1. Ya vimos los resultados que da.

  2. En la lista de $iNextEst se dan dos resultados:

  1. $size_based (LCL y UCL se obtienen para tamaño fijo). Para el tamaño de muestra correspondiente a cada nudo, el primer marco de datos (como se muestra en $size_based) incluye el nombre del Ensamblaje, el tamaño de la muestra (m, es decir, cada uno de los 40 nudos), el método (Rarefacción, Observado o Extrapolación). , dependiendo de si el tamaño m es menor, igual o mayor que el tamaño de muestra de referencia), el orden de diversidad (orden.q), la estimación de diversidad de orden q (qD), los límites de confianza superior e inferior del 95% de diversidad (qD.LCL, qD.UCL) y la estimación de cobertura de la muestra (SC) junto con los límites de confianza superior e inferior del 95% de la cobertura de la muestra (SC.LCL, SC.UCL). Estas estimaciones de cobertura muestral con intervalos de confianza se utilizan para trazar la curva de integridad de la muestra.
  2. $coverage_based (LCL y UCL se obtienen para cobertura fija; la longitud del intervalo es más amplia debido a la variación del tamaño en los bootstraps). Da los mismos resultados que el $size_based solo que variarán/diferirán dependiendo de la cantidad de bootstraps que se pongan.
  1. $AsyEst enlista el nombre del “Ensamblaje”, la diversidad (riqueza de especies para q = 0, diversidad de Shannon para q = 1 y diversidad de Simpson para q = 2), la diversidad observada (Observed), la estimación de diversidad asintótica (Estimator), el “s.e.” del estimador asintótico (s.e.) y los límites de confianza superior e inferior del 95% asociados (LCL, UCL). Los estimadores de diversidad como Shannon y Simpson se transforman con este comando.

Para graficar los resultados estadísticos obtenidos en iNEXT

ggiNEXT(x, type=1, se=TRUE, facet.var=“None”, color.var=“Assemblage”, grey=FALSE)

Los “type” son dependiendo de los resultados que necesitemos: 1 = para diversidad vs número de individuos, 2 = cobertura vs número de individuos, 3 = diversidad vs cobertura. Con “se” (intervalos de confianza) por defecto es = True. El “facet.var” nos dice como agrupar los plots, si hacer solo 1, separarlo por ensamblaje, etc. El “color.var” nos dice qué color deseamos para la curva de resultados, puede variar entre los sitios, pero mejor dejarlas del mismo color para ambos gráficos El “grey” = FALSE evita que se ponga el gráfico en escala de grises, está por defecto.

Ejemplo

Primero llamemos a los paquetes que necesitariamos

library(iNEXT)
library(ggplot2)
library(gridExtra)
library(grid)

Utilizaremos la base de datos “bird” que contiene la abundancia de especies de aves en dos sitios (Norte y Sur) del Barrington Tops National Park, Australia (Chao et al. 2015).

data(bird)
str(bird)
## 'data.frame':    41 obs. of  2 variables:
##  $ North.site: int  0 0 41 0 3 1 5 4 4 11 ...
##  $ South.site: int  3 18 31 2 1 2 5 1 6 32 ...
ausbird <- iNEXT(bird, q=c(0, 1, 2), datatype="abundance", endpoint = 500)

Diversidad vs número de individuos:

ausbird.ggp1 <- ggiNEXT(ausbird, type=1, facet.var="Assemblage", color.var = "Assemblage")
ausbird.ggp1 + theme_classic() + theme(legend.position = "bottom", legend.box = "horizontal")

En el eje X: número de individuos, en el eje Y: diversidad de especies. anaranjado (0) indica riqueza de especies, azul (1) indica índice de diversidad de Shannon y rosado (2) indica índice de Simpson. (0): Para el sitio norte, la mayor riqueza de especies se alcanzó a los 200 individuos muestreados, mostrando una diversidad de 28 especies. Para el sitio sur, la mayor diversidad de alcanzó a los 300 individuos muestreados, mostrando una diversidad de 40 especies. (1): El índice de Shannon indica que mientras el valor sea más alto, mayor será la diversidad de un sitio. En este caso, el sitio más diverso fue el sur. (2): El índice de Simpson indica que mientras más pequeño sea un valor, más diverso será el sitio. En este caso, el sitio norte muestra un índice menor al sitio sur.

Para comprobar que los datos del plot coinciden con nuestros índices, corremos los comandos que vimos al inicio pero con la transformación:

ChaoRichness(bird, datatype = "abundance", conf = 0.95)
##            Observed Estimator Est_s.e. 95% Lower 95% Upper
## North.site       27    31.478    4.781    27.810    51.746
## South.site       38    40.077    2.499    38.325    51.267
ChaoShannon(bird, datatype = "abundance", transform = TRUE, conf = 0.95, B = 200)
##            Observed Estimator Est_s.e 95% Lower 95% Upper
## North.site   16.536    17.966   1.265    16.536    20.446
## South.site   25.447    27.251   1.303    25.447    29.805
ChaoSimpson(bird, datatype = "abundance", transform = TRUE, conf = 0.95, B = 200)
##            Observed Estimator Est_s.e. 95% Lower 95% Upper
## North.site   11.848    12.524    1.385    11.848    15.238
## South.site   19.639    20.913    1.407    19.639    23.671

Cobertura vs número de individuos:

ausbird.ggp2 <- ggiNEXT(ausbird, type=2, facet.var="Assemblage", color.var = "Assemblage")
ausbird.ggp2 + theme_classic() + theme(legend.position = "bottom", legend.box = "horizontal")

Eje X: número de individuos, eje Y: completitud de muestra. El sample coverage o completitud de la muestra se interpreta como la cantidad de individuos presentes en un sitio dividido la cantidad de individuos observados en el sitio. Entre más completa o (mejor cubierta) se tenga la muestra, es menos probable que un individuo colectado en el próximo muestreo sea de una especie no resgistrada anteriormente.

En este caso, para ambos sitios existe una completitud de muestra similar, ya que ambos arrojan un valor de 1.

Diversidad vs. cobertura:

ausbird.ggp3 <- ggiNEXT(ausbird, type=3, facet.var="Assemblage", color.var = "Assemblage")
ausbird.ggp3 + theme_classic() + theme(legend.position = "bottom", legend.box = "horizontal")

Eje X: completitud de la muestra, eje Y: diversidad de especies. Este gráfico representa el inverso del primer gráfico. Para mabos sitios es posible observar que la diversidad de especies está directamente relacionada a la completitud de la muestra, ya que la representatividad de una muestra depende no solo del número de especies que faltan por registrar, sino de sus abundancias promedio.

Glosario

Rarefacción:

Método que se utiliza para comparar índices de diversidad entre sitios en base a un mismo número de individuos. El sitio con mayor numero de individuos se submuestrea sin reemplazo aleatoriamente con múltiples ejecuciones para generar un índice promedio que se puede comparar con el índice de otro sitio en base a un mismo número de individuos (Carmona-Galindo y Carmona, 2013).

Interpolación:

Método en el que se utilizan valores conocidos para estimar un valor desconocido (Carmona-Galindo y Carmona, 2013).

Extrapolación:

Proceso de estimar más allá del intervalo de observación original el valor de la variable,con base en su relación con otra variable. Está sujeta a un grado alto de incertidumbre y riesgo de producir resultados no significativos (Carmona-Galindo y Carmona, 2013).

Abundancia:

Representación relativa de una especie en un sitio específico.Se mide como el número de individuos encontrados en una muestra (Chao et al. 2014).

Incidencia:

Número de eventos en el cual aparece una especie en un período determinado (Chao et al. 2014).

Hill numbers:

También llamado número efectivo de especies. Se utilizan para caracterizar la diversidad taxonómica, filogenética o funcional de una comunidad o sitio. Las estimaciones empíricas de los Hill numbers (incluída la riqueza de especies) tienden a ser una función del esfuerzo de muestreo y, por lo tanto, tienden a aumentar a medida que la muestra se completa (Chao et al. 2014).

Referencias

Carmona-Galindo, V.D., & Carmona, T.V. (2013). La Diversidad de los Analisis de Diversidad. Bioma, 14, 20-28.

Chao, A. (08 de noviembre de 2023). iNEXT. Anne Chao’s Website. http://chao.stat.nthu.edu.tw/wordpress/software_download/inext-online/

Chao, A., Gotelli, N. J., Hsieh, T. C., Sander, E. L., Ma, K. H., Colwell, R. K., & Ellison, A. M. (2014). Rarefaction and extrapolation with Hill numbers: a framework for sampling and estimation in species diversity studies. Ecological monographs, 84(1), 45-67.

Escalante, T. (2003). ¿ Cuántas especies hay? Los estimadores no paramétricos de Chao. Elementos 52: 53–56.

Hsieh, T.C., Ma, K.H. & Chao, A. (2016) iNEXT: An R package for interpolation and extrapolation of species diversity (Hill numbers). Methods in Ecology and Evolution, 7, 1451-1456.

Manual disponible en Rpubs: https://rpubs.com/AnaluArevaloF/1112645