Este documento describe el análisis de la vegetación en dos tipos de hábitats imaginarios, que llamé cerca y lejos, indicando si están cerca o lejos del camino. Los datos se ‘tomaron’ en parcelas o transectas de igual tamaño. Se ‘midieron’ tres variables ambientales en cada parcela: porcentaje de terreno rocoso (rocoso), porcentaje del sotobosque con plantas invasoras (invasion), pendiente en grados (pendiente). Utizo el paquete picante para los análisis.
library(picante)
library(gt)
Los datos están en dos archivos .csv (comma separated values) producidos a partir de hojas de cálculo Excel. Un archivo es de abundancia en una matriz de parcelas X especies, y la otra es una tabla con las medidas de las variables ambientales por parcela.
# abundancia
comm <- read.csv("diversidad.comunidades.imagina.csv",
header = TRUE, row.names = 1)
head(comm)
## sp_1 sp_2 sp_3 sp_4 sp_5 sp_6 sp_7 sp_8 sp_9 sp_10 sp_11 sp_12 sp_13
## A_1 0 3 1 10 4 0 0 1 6 1 1 1 0
## B_1 1 1 0 2 1 0 1 0 0 0 1 0 0
## A_2 1 2 1 3 0 1 0 2 3 1 0 1 0
## B_2 2 1 1 4 2 0 1 1 2 0 1 1 0
## A_3 0 3 1 3 1 0 2 1 5 2 1 3 0
## B_3 3 1 0 1 2 0 3 0 1 1 0 0 3
# ambientales
metadata <- read.csv("plot.metadata.imagina.csv",
header = TRUE, row.names = 1)
head(metadata)
## elevacion pendiente rocoso invasion
## A_1 7 20 50 75
## B_1 7 0 25 100
## A_2 8 0 0 50
## B_2 8 10 25 25
## A_3 10 10 50 25
## B_3 10 20 75 25
En esta sección mostramos el análisis de la diversidad.
Gráficas de acumulación de especies, en el estudio.
# conjunta
plot(specaccum(comm), xlab = "# de muestras", ylab = "# de especies")
A continuación calculamos índices de diversidad para cada hábitat.
# índices de biodiversidad: H:Shannon-Wiener / S:# especies / J:igualdad Pielou
# area 1
area1 <- comm[1:2, ]
sum.area1 <- apply(area1, 2, sum)
Ha1 <- diversity(sum.area1, index = "shannon", MARGIN = 1, base = exp(1))
Sa1 <- specnumber(sum.area1)
Ja1 <- Ha1/log(Sa1)
indi.a1 <- c(Sa1,Ha1,Ja1)
# area 2
area2 <- comm[3:4, ]
sum.area2 <- apply(area2, 2, sum)
Ha2 <- diversity(sum.area2, index = "shannon", MARGIN = 1, base = exp(1))
Sa2 <- specnumber(sum.area2)
Ja2 <- Ha2/log(Sa2)
indi.a2 <- c(Sa2,Ha2,Ja2)
# area 3
area3 <- comm[5:6, ]
sum.area3 <- apply(area3, 2, sum)
Ha3 <- diversity(sum.area3, index = "shannon", MARGIN = 1, base = exp(1))
Sa3 <- specnumber(sum.area3)
Ja3 <- Ha3/log(Sa3)
indi.a3 <- c(Sa3,Ha3,Ja3)
# area 4
area4 <- comm[7:8, ]
sum.area4 <- apply(area4, 2, sum)
Ha4 <- diversity(sum.area4, index = "shannon", MARGIN = 1, base = exp(1))
Sa4 <- specnumber(sum.area4)
Ja4 <- Ha2/log(Sa4)
indi.a4 <- c(Sa4,Ha4,Ja4)
# tabla
indices <- rbind(indi.a1, indi.a2, indi.a3, indi.a4)
rownames(indices) <- c("Area 1", "Area 2", "Area 3", "Area 4")
indices.gt <- as.data.frame(indices)
indices.gt %>%
gt(rownames_to_stub = TRUE) %>%
tab_header(
title = "Tabla 1. Indices de biodiversidad para las áreas"
) %>%
cols_label(
V1 = html("Riqueza, S"),
V2 = html("Shannon, H"),
V3 = html("Equidad, J")
) %>%
fmt_number(
columns = 3:4,
decimals = 2
)
Tabla 1. Indices de biodiversidad para las áreas | |||
---|---|---|---|
Riqueza, S | Shannon, H | Equidad, J | |
Area 1 | 11 | 1.97 | 0.82 |
Area 2 | 12 | 2.28 | 0.92 |
Area 3 | 12 | 2.36 | 0.95 |
Area 4 | 11 | 2.13 | 0.95 |
Vamos a construir un agrupamiento jerárquico de acuerdo a la composición de especies, basado en las disimilaridades Bray-Curtis.
# hierachircal clustering - agrupamiento jerárquico
# cálculo de distancias Bray-Curtis entre muestras
comm.bc.dist <- vegdist(comm, 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 - disimilaridades")
Estaremos usando el procedimiento Non Metric MultiDimensional. NMDS es una técnica utilizada para simplificar los datos multivariados en unos pocos ejes para facilitar el reconocimiento e interpretación de patrones y diferencias entre grupos.
Ordenación de las parcelas (plots)
# La función metaMDS transforma automáticamente los datos y verifica la robustez
# de la solución
# robustez
comm.bc.mds <- metaMDS(comm, dist = "bray")
## Wisconsin double standardization
## Run 0 stress 0.04086855
## Run 1 stress 0.1658918
## Run 2 stress 0.04086855
## ... Procrustes: rmse 5.151469e-06 max resid 8.452147e-06
## ... Similar to previous best
## Run 3 stress 0.1389373
## Run 4 stress 0.04086856
## ... Procrustes: rmse 4.269765e-05 max resid 6.796546e-05
## ... Similar to previous best
## Run 5 stress 0.1389373
## Run 6 stress 0.04086856
## ... Procrustes: rmse 2.030173e-05 max resid 3.230912e-05
## ... Similar to previous best
## Run 7 stress 0.1358949
## Run 8 stress 0.1658918
## Run 9 stress 0.04086856
## ... Procrustes: rmse 2.158493e-05 max resid 3.433225e-05
## ... Similar to previous best
## Run 10 stress 0.1389375
## Run 11 stress 0.1389374
## Run 12 stress 0.2852161
## Run 13 stress 0.04086856
## ... Procrustes: rmse 3.71363e-05 max resid 5.87995e-05
## ... Similar to previous best
## Run 14 stress 0.04086858
## ... Procrustes: rmse 6.091492e-05 max resid 9.696955e-05
## ... Similar to previous best
## Run 15 stress 0.1658918
## Run 16 stress 0.04086855
## ... Procrustes: rmse 7.795419e-06 max resid 1.038065e-05
## ... Similar to previous best
## Run 17 stress 0.04086858
## ... Procrustes: rmse 4.56656e-05 max resid 7.188633e-05
## ... Similar to previous best
## Run 18 stress 0.2057087
## Run 19 stress 0.04086855
## ... New best solution
## ... Procrustes: rmse 2.803409e-05 max resid 4.593094e-05
## ... Similar to previous best
## Run 20 stress 0.04086855
## ... New best solution
## ... Procrustes: rmse 2.210517e-05 max resid 3.639587e-05
## ... Similar to previous best
## *** Solution reached
# plot site scores as text
ordiplot(comm.bc.mds, display = "sites", type = "text")
Incluyendo las especies
# con especies
# automated plotting of results - tries to eliminate overlapping labels
ordipointlabel(comm.bc.mds, cex = (c(1, 1)))
Grafica simplificada de parcelas
# plot anything yet
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 = "blue", select = metadata$elevacion ==
7)
points(mds.fig, "sites", pch = 19, col = "green", select = metadata$elevacion ==
8)
points(mds.fig, "sites", pch = 19, col = "orange", select = metadata$elevacion ==
10)
points(mds.fig, "sites", pch = 19, col = "red", select = metadata$elevacion ==
14)
Variables ambientales
# parametros ambientales
# ordiplot(comm.bc.mds)
# calculate and plot environmental variable correlations with the axes use
# the subset of metadata that are environmental data
ordipointlabel(comm.bc.mds, cex = (c(1, 1)))
par
## function (..., no.readonly = FALSE)
## {
## .Pars.readonly <- c("cin", "cra", "csi", "cxy", "din", "page")
## single <- FALSE
## args <- list(...)
## if (!length(args))
## args <- as.list(if (no.readonly)
## .Pars[-match(.Pars.readonly, .Pars)]
## else .Pars)
## else {
## if (all(unlist(lapply(args, is.character))))
## args <- as.list(unlist(args))
## if (length(args) == 1) {
## if (is.list(args[[1L]]) || is.null(args[[1L]]))
## args <- args[[1L]]
## else if (is.null(names(args)))
## single <- TRUE
## }
## }
## value <- .External2(C_par, args)
## if (single)
## value <- value[[1L]]
## if (!is.null(names(args)))
## invisible(value)
## else value
## }
## <bytecode: 0x55e09d1189d8>
## <environment: namespace:graphics>
plot(envfit(comm.bc.mds, metadata[, 1:4]))