Introducción

Para este ejercicio vamos a analizar los datos del proyecto de biodiversidad en la playa norte (bahía de Seven Seas) de la Reserva Cabezas de San Juan, Fajardo: Datos de Biodiversidad de la Playa.
Las muestras (20) se tomaron a lo largo de la zona de eoleanita de la playa y de la pradera de Thalassia tertudinum. Las unidades de muestreos fueron cuadrados de 0.5 m (0.25 \(m^2\)). Todos los organismos fueron contados (abundancia) o estimada su cobertura (en porciento).

Cargar los datos de la playa

A continuación un ejemplo de lectura de los datos del archivo en Excel. Cargamos primero los datos de la zona de eoleanita (Eole) y pradera de Thalassia (Thal) por separado y luego los datos de ambas comunidades.

# paquete para leer directamente datos en Excel
library(readxl)
# paquete para 'limpiar' los datos (nombres de las columnas/especies)
library(janitor)
eoleanita <- read_excel("data/LabEcoCom biodiversidad playa 2023.xlsx", 
                       sheet = "eoleanita")
# limpiando el data frame (nombres de especies)
eoleanita <- clean_names(eoleanita) # ver contenido en Global Environment

Cambiar datos NA por 0

En los datos originales tenemos espacios en blanco (NA) para cuando no se encontraron individuos de una especie; esto puede causar problemas en los análisis, así que vamos a cambiar esos NA por ceros (0).

# cambio de NA a 0
eoleanita <- data.frame(lapply(eoleanita, function(x) ifelse(is.na(x), 0, x)))

Paquetes necesarios para el análisis de biodiversidad

BiodiversityR puede presentar problemas…

library(vegan)
library(BiodiversityR)
library(gt)

Indices de Biodiversidad por muestra

Calcular los índices de Shannon-Wiener (H), Simpson (D), riqueza de especies (S), y equidad (‘evenness’, J), en cada bráctea, y para cada año.

#usando vegan
H <- diversity(eoleanita, index = "shannon")
D <- diversity(eoleanita, index = "simpson")
S <- specnumber(eoleanita)
J <- H/log(specnumber(eoleanita))
indices <- data.frame(S,H,D,J)
gt(indices)
S H D J
4 1.04506884 0.59760000 0.7538578
3 0.80090270 0.46713240 0.7290131
5 1.34854893 0.70291637 0.8379006
4 1.09855938 0.64882414 0.7924431
3 0.90879543 0.56703398 0.8272213
3 0.56467367 0.30399700 0.5139881
4 0.88985488 0.54518430 0.6418946
4 1.23507581 0.67969136 0.8909189
4 0.83998578 0.52531200 0.6059217
3 0.60016651 0.32622268 0.5462951
3 0.36764948 0.16942149 0.3346490
4 0.21655989 0.07993995 0.1562149
4 0.98192876 0.57140362 0.7083119
4 0.78682158 0.40875000 0.5675718
5 1.36868964 0.70312500 0.8504147
3 0.78670124 0.47530864 0.7160863
2 0.07851576 0.02984389 0.1132743
4 0.54586925 0.25632653 0.3937614
3 0.83052369 0.50781250 0.7559752
1 0.00000000 0.00000000 NaN

Estadísticas básicas de los índices

A partir de los resultados anteriores (índices de diversidad y equidad, para cada muestra) puede obtener algunas estadísticas descriptivas, incluyendo parámetros estadísticos como media, mediana, desviación estándar, et c., para cada uno de ellos.

med <- mean(H)
var <- var(H)
max <- max(H)
min <- min(H)
indxstatH <- data.frame(med,var,max,min)
gt(indxstatH)
med var max min
0.7647446 0.1495577 1.36869 0

Indices en conjunto (“pooled”)

Además de calcular los índices de diversidad para cada cuadrado, puedes calcular los índices tomando todos los datos en conjunto:

#índices de diversidad en conjunto con BiodiversityR
Sp <- diversityresult(eoleanita, index=("richness"), method=("pooled"))
Hp <- diversityresult(eoleanita, index=("Shannon"), method=("pooled"))
Dp <- diversityresult(eoleanita, index=("Simpson"), method=("pooled"))
Jp <- diversityresult(eoleanita, index=("Jevenness"), method=("pooled"))
indxpool <- data.frame(Sp[1,1],Hp[1,1],Dp[1,1],Jp[1,1])
gt(indxpool)
Sp.1..1. Hp.1..1. Dp.1..1. Jp.1..1.
23 2.2311136 0.85873512 0.7115668

Curvas de acumulación de especies

Podemos estudiar la acumulación de especies (nuevas) en relación al número de muestras, o al número de individuos seleccionados, en ambos casos, mediante un proceso ‘aleatorio’.

sac <- specaccum(eoleanita, method = "collector")
plot(sac, ci.type="polygon", ci.col="yellow") #ver vegan para opciones

sac <- specaccum(eoleanita, method = "exact")
plot(sac, ci.type="polygon", ci.col="green")

#por individuos
sac <- specaccum(eoleanita, method = "rarefaction")
plot(sac, xvar = "individual", ci.type="polygon", ci.col="yellow") #ver vegan para opciones

Rarefacción

El método de rarefacción, además de servir para construir las curvas especies-individuos, permite estandarizar las comparaciones de las curvas de especie-individuos, y estimar la riqueza de especies, ajustando las muestras al valor mínimo de individuos encontrados en las muestras (o uno asignado).

# curvas especies - muestra
rarecurve(eoleanita, step = 1, col = "blue", cex = 0.6)

Curvas abundancia vs rango

Otro aspecto importante en el análisis de la biodiversidad, es determinar como se distribuyen las especies, de acuerdo a su abundancia. Así podemos determinar especies dominantes, subdominantes, en crecimiento, especies raras, et c.

eoleanitadf <- as.data.frame(eoleanita) #para forzar al formato data frame
# usando BiodiversityR
RkAb <- rankabundance(eoleanitadf)
head(RkAb,10)
##             rank abundance proportion plower pupper accumfreq logabun rankfreq
## algaverde_2    1       408       22.6    8.2   37.0      22.6     2.6      3.0
## alga_1         2       380       21.0    8.7   33.4      43.6     2.6      6.1
## caracol_3      3       280       15.5    1.2   29.8      59.1     2.4      9.1
## poliqueto_1    4       138        7.6   -1.2   16.4      66.7     2.1     12.1
## alga_3         5       130        7.2   -3.2   17.6      73.9     2.1     15.2
## algaroja_2     6       128        7.1   -2.0   16.2      81.0     2.1     18.2
## padina_1       7        88        4.9   -2.9   12.7      85.9     1.9     21.2
## algaverde_1    8        81        4.5   -1.2   10.1      90.4     1.9     24.2
## alga_2         9        47        2.6   -2.0    7.2      93.0     1.7     27.3
## caracol_1     10        43        2.4   -0.4    5.2      95.4     1.6     30.3
rankabunplot(RkAb, scale='proportion', addit=FALSE, specnames=c(1,2,3))

Análisis de Biodiversidad entre Comunidades

Vamos a utilizar algunos de los procedimientos presentados por Kembel (2013).

# paquete con sus dependencias
library(picante)
# para manejar los datos
library(dplyr)
# para gráficas
library(ggplot2)

Datos

Los datos deben estar en el siguiente formato:
primera fila

muestra especie_1 especie_2 … especie_n
m1 c1 c2 … cn
…

donde ci es la cantidad de individuos de la especie i para cada muestra.

Entrada de datos al sistema

comunidades <- read.csv("data/LabEcoCom biodiversidad playa 2023.csv", header = TRUE, row.names = 1)
# cambio de NA a 0
comunidades[is.na(comunidades)] <- 0
# revisar si es un data.frame y sus dimensiones
str(comunidades)
## 'data.frame':    39 obs. of  33 variables:
##  $ caracol_1             : num  19 9 0 0 0 0 0 0 0 9 ...
##  $ caracol_2             : num  2 0 0 1 0 0 0 0 0 0 ...
##  $ algaroja_1            : num  4 0 0 0 0 0 0 0 0 0 ...
##  $ caracol_3             : num  0 25 151 44 0 0 0 0 0 0 ...
##  $ alga_1                : num  25 75 80 70 0 0 0 80 0 50 ...
##  $ algaverde_1           : num  0 0 10 0 0 0 0 0 50 3 ...
##  $ alga_2                : num  0 0 0 40 0 0 0 0 0 0 ...
##  $ padina_1              : num  0 0 80 0 0 0 2 0 1 0 ...
##  $ algaverde_2           : num  0 0 30 0 30 3 0 40 0 0 ...
##  $ poliqueto_1           : num  0 0 0 0 31 0 0 47 0 0 ...
##  $ algaverde_3           : num  0 0 0 0 0 0 4 0 0 0 ...
##  $ algaroja_2            : num  0 0 0 0 5 60 50 13 0 0 ...
##  $ algaverde_4           : num  0 0 0 0 0 10 0 0 0 0 ...
##  $ alga_3                : num  0 0 0 0 0 0 60 0 70 0 ...
##  $ Valonia_aegagropila   : num  0 0 0 0 0 0 0 0 4 0 ...
##  $ algaverde_6           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ algaverde_7           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ algaverde_8           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ nudibranquio_1        : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ equinodermo_1         : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ algaverde_9           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Hermodice_carunculata : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ caracol_4             : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Penicillus_dumetosus  : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ coral_1               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ esponja_1             : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Penicillus_pyriformis : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ algaverde_10          : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ algaparda_1           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ bivalvo_1             : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ alga_5                : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ alga_6                : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Cleisostoma_filiformis: num  0 0 0 0 0 0 0 0 0 0 ...

Revisamos el total de individuos por muestra (abundancia), para verificar que no haya muestras con cero datos. También estandarizamos los datos a abundancia relativa en cada muestra, usando la función decostand del paquete vegan.

library(vegan)
# verificar abundancia por muestra
apply(comunidades, 1, sum)
##  Eole_1  Eole_2  Eole_3  Eole_4  Eole_5  Eole_6  Eole_7  Eole_8  Eole_9 Eole_10 
##      50     109     351     155      66      73     116     180     125      62 
## Eole_11 Eole_12 Eole_13 Eole_14 Eole_15 Eole_16 Eole_17 Eole_18 Eole_19 Eole_20 
##      44      73     107      40      16      90      66      35      48       1 
##  Thal_2  Thal_3  Thal_4  Thal_5  Thal_6  Thal_7  Thal_8  Thal_9 Thal_10 Thal_11 
##       2       1       1       1       8       5       2       8       9       4 
## Thal_12 Thal_13 Thal_14 Thal_15 Thal_16 Thal_17 Thal_18 Thal_19 Thal_20 
##      47      23      10       1       3      12      23      17      29
# estandarizar los datos por muestra
comunidades <- decostand(comunidades, method = "total")
# verificamos estandarización
apply(comunidades, 1, sum)
##  Eole_1  Eole_2  Eole_3  Eole_4  Eole_5  Eole_6  Eole_7  Eole_8  Eole_9 Eole_10 
##       1       1       1       1       1       1       1       1       1       1 
## Eole_11 Eole_12 Eole_13 Eole_14 Eole_15 Eole_16 Eole_17 Eole_18 Eole_19 Eole_20 
##       1       1       1       1       1       1       1       1       1       1 
##  Thal_2  Thal_3  Thal_4  Thal_5  Thal_6  Thal_7  Thal_8  Thal_9 Thal_10 Thal_11 
##       1       1       1       1       1       1       1       1       1       1 
## Thal_12 Thal_13 Thal_14 Thal_15 Thal_16 Thal_17 Thal_18 Thal_19 Thal_20 
##       1       1       1       1       1       1       1       1       1

Metadata

Creamos los nombres de las comunidades, para poder luego realizar comparaciones.

metadata <- read.csv("data/comunidades-playa.csv", header = TRUE, row.names = 1)
metadata
##         comunidades
## Eole_1    eoleanita
## Eole_2    eoleanita
## Eole_3    eoleanita
## Eole_4    eoleanita
## Eole_5    eoleanita
## Eole_6    eoleanita
## Eole_7    eoleanita
## Eole_8    eoleanita
## Eole_9    eoleanita
## Eole_10   eoleanita
## Eole_11   eoleanita
## Eole_12   eoleanita
## Eole_13   eoleanita
## Eole_14   eoleanita
## Eole_15   eoleanita
## Eole_16   eoleanita
## Eole_17   eoleanita
## Eole_18   eoleanita
## Eole_19   eoleanita
## Eole_20   eoleanita
## Thal_2    thalassia
## Thal_3    thalassia
## Thal_4    thalassia
## Thal_5    thalassia
## Thal_6    thalassia
## Thal_7    thalassia
## Thal_8    thalassia
## Thal_9    thalassia
## Thal_10   thalassia
## Thal_11   thalassia
## Thal_12   thalassia
## Thal_13   thalassia
## Thal_14   thalassia
## Thal_15   thalassia
## Thal_16   thalassia
## Thal_17   thalassia
## Thal_18   thalassia
## Thal_19   thalassia
## Thal_20   thalassia
# verificamos nombres de filas en comunidades y metadata
all.equal(rownames(comunidades), rownames(metadata))
## [1] TRUE

Análisis multivariado para similaridad entre comunidades

Dendograma de grupos

Clasificación de las muestras en grupos (‘clusters’) de acuerdo a su similitud de composición de especies.

# calculate Bray-Curtis distance among samples
comm.bc.dist <- vegdist(comunidades, method = "bray")
# cluster communities using average-linkage algorithm
comm.bc.clust <- hclust(comm.bc.dist, method = "average")
# plot cluster diagram
plot(comm.bc.clust, ylab = "Bray-Curtis dissimilarity")

Ordenación

Arreglo de las muestras y taxones en dos ejes para máxima separación o agrupación.

# The metaMDS function automatically transforms data and checks solution
# robustness
comm.bc.mds <- metaMDS(comunidades, dist = "bray")
## Run 0 stress 0.066974921 
## Run 1 stress 0.074301877 
## Run 2 stress 0.073297741 
## Run 3 stress 0.067460171 
## ... Procrustes: rmse 0.069479339  max resid 0.19505114 
## Run 4 stress 0.074098791 
## Run 5 stress 0.066514035 
## ... New best solution
## ... Procrustes: rmse 0.03378727  max resid 0.1523122 
## Run 6 stress 0.071408688 
## Run 7 stress 0.064215163 
## ... New best solution
## ... Procrustes: rmse 0.056461718  max resid 0.22435365 
## Run 8 stress 0.065196953 
## Run 9 stress 0.066981629 
## Run 10 stress 0.071408765 
## Run 11 stress 0.071819406 
## Run 12 stress 0.070589905 
## Run 13 stress 0.071269362 
## Run 14 stress 0.078350527 
## Run 15 stress 0.07140833 
## Run 16 stress 0.067592977 
## Run 17 stress 0.084262567 
## Run 18 stress 0.083667746 
## Run 19 stress 0.06639722 
## Run 20 stress 0.071168116 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      4: no. of iterations >= maxit
##     16: stress ratio > sratmax
# Assess goodness of ordination fit (stress plot)
stressplot(comm.bc.mds)

# plot site scores as text
ordiplot(comm.bc.mds, display = "sites", type = "text")

# automated plotting of results - tries to eliminate overlapping labels
ordipointlabel(comm.bc.mds)
# ordination plots are highly customizable set up the plotting area but
# don't plot anything yet
title("Figura O1")

mds.fig <- ordiplot(comm.bc.mds, type = "none")
# plot just the samples, colour by habitat, pch=19 means plot a circle
points(mds.fig, "sites", pch = 19, col = "green", 
       select = metadata$comunidad == "bosque")
points(mds.fig, "sites", pch = 19, col = "blue", 
       select = metadata$comunidad == "plantacion")
# add confidence ellipses around habitat types
title("Figura O2")
ordiellipse(comm.bc.mds, metadata$comunidad, conf = 0.95, label = TRUE)
# overlay the cluster results we calculated earlier
ordicluster(comm.bc.mds, comm.bc.clust, col = "gray")

Gradiente de abundancia de los taxones

Gráfica de isolíneas de fracción de abundancia de los taxones.

# plot alga_3 abundance. cex increases the size of bubbles.
ordisurf(comm.bc.mds, comunidades[, "alga_3"], bubble = TRUE, main = "Abundancia de alga_3", 
    cex = 3)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
## 
## Estimated degrees of freedom:
## 7.12  total = 8.12 
## 
## REML score: -8.9738442
# plot alga_1 abundance. cex increases the size of bubbles.
ordisurf(comm.bc.mds, comunidades[, "alga_1"], bubble = TRUE, main = "Abundancia de alga_1", 
    cex = 3)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
## 
## Estimated degrees of freedom:
## 5.86  total = 6.86 
## 
## REML score: -14.028413