1 Introducción a la Diversidad Beta

1.1 ¿Qué es la diversidad Beta?

La diversidad beta (β-diversidad) es un concepto fundamental en ecología de comunidades que cuantifica la variación en la composición de especies entre diferentes sitios, hábitats o comunidades. Fue introducida por Robert Whittaker en 1960 como parte de su esquema jerárquico de diversidad, que incluye:

  • Diversidad alfa (α): Diversidad local dentro de un sitio
  • Diversidad beta (β): Diversidad entre sitios (recambio de especies)
  • Diversidad gamma (γ): Diversidad regional total

1.2 Importancia ecológica

La diversidad beta es crucial para entender:

  • Procesos de especiación y extinción
  • Ensamblaje de comunidades
  • Efectos de la fragmentación del hábitat
  • Patrones biogeográficos
  • Conservación de la biodiversidad

2 Fundamentos Teóricos

2.1 Conceptos básicos de Whittaker

Whittaker definió la diversidad beta como la relación entre la diversidad gamma (regional) y la diversidad alfa (local):

\[\beta = \frac{\gamma}{\alpha}\]

Esta formulación multiplicativa indica cuántas veces mayor es la diversidad regional comparada con la diversidad local promedio.

2.2 Formulación moderna

En el contexto moderno, la diversidad beta se entiende mejor como una medida de diferenciación o disimilitud entre comunidades. Existen dos enfoques principales:

2.2.1 1. Enfoque de disimilitud

Basado en índices que van de 0 (comunidades idénticas) a 1 (comunidades completamente diferentes).

2.2.2 2. Enfoque multiplicativo

Mantiene la interpretación original de Whittaker pero con refinamientos metodológicos.

2.3 Componentes de la diversidad Beta

La diversidad beta puede descomponerse en dos componentes principales:

  • Recambio (turnover): Reemplazo de especies entre sitios
  • Anidamiento (nestedness): Diferencias en riqueza donde sitios pobres son subconjuntos de sitios ricos

3 Índices y Métricas de Diversidad Beta

3.1 Índices clásicos de disimilitud

3.1.1 Índice de Sørensen

\[\beta_{sor} = \frac{b + c}{2a + b + c}\]

Donde: - a = número de especies compartidas - b = especies únicas del sitio 1 - c = especies únicas del sitio 2

3.1.2 Índice de Jaccard

\[\beta_{jac} = \frac{b + c}{a + b + c}\]

3.1.3 Índice de Bray-Curtis

Para datos de abundancia:

\[BC = \frac{\sum_{i=1}^{n}|x_{i} - y_{i}|}{\sum_{i=1}^{n}(x_{i} + y_{i})}\]

3.2 Enfoque de Baselga

Baselga (2010) propuso una descomposición aditiva de la diversidad beta:

\[\beta_{sor} = \beta_{sim} + \beta_{nes}\]

Donde: - βₛᵢₘ = componente de recambio (Simpson) - βₙₑₛ = componente de anidamiento

# Cargar librerías necesarias
library(vegan)
## Loading required package: permute
## Loading required package: lattice
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(betapart)
library(adespatial)
## Registered S3 methods overwritten by 'adegraphics':
##   method         from
##   biplot.dudi    ade4
##   kplot.foucart  ade4
##   kplot.mcoa     ade4
##   kplot.mfa      ade4
##   kplot.pta      ade4
##   kplot.sepan    ade4
##   kplot.statis   ade4
##   scatter.coa    ade4
##   scatter.dudi   ade4
##   scatter.nipals ade4
##   scatter.pco    ade4
##   score.acm      ade4
##   score.mix      ade4
##   score.pca      ade4
##   screeplot.dudi ade4
## Registered S3 method overwritten by 'spdep':
##   method   from
##   plot.mst ape
## Registered S3 methods overwritten by 'adespatial':
##   method             from       
##   plot.multispati    adegraphics
##   print.multispati   ade4       
##   summary.multispati ade4
library(ade4)
## 
## Attaching package: 'ade4'
## The following object is masked from 'package:adespatial':
## 
##     multispati
library(corrplot)
## corrplot 0.92 loaded
library(knitr)
library(DT)

4 Ejemplo conceptual de descomposición de Baselga

5 Crear datos sintéticos

set.seed(123) site1 <- c(1, 1, 1, 0, 0, 0, 0) site2 <- c(1, 1, 0, 1, 1, 0, 0) site3 <- c(1, 0, 0, 0, 0, 1, 1)

6 Crear matriz de comunidades

comm_matrix <- rbind(site1, site2, site3) colnames(comm_matrix) <- paste(“sp”, 1:7, sep=““) rownames(comm_matrix) <- c(”Sitio 1”, “Sitio 2”, “Sitio 3”)

print(“Matriz de comunidades de ejemplo:”) kable(comm_matrix)

7 Visualizar matriz de comunidades

library(reshape2) comm_melted <- melt(comm_matrix) colnames(comm_melted) <- c(“Sitio”, “Especie”, “Presencia”)

ggplot(comm_melted, aes(x = Especie, y = Sitio, fill = factor(Presencia))) + geom_tile(color = “white”, size = 0.5) + scale_fill_manual(values = c(“0” = “white”, “1” = “darkblue”)) + theme_minimal() + labs(title = “Matriz de presencia-ausencia de especies”, fill = “Presencia”) + theme(axis.text.x = element_text(angle = 45, hjust = 1))


# Protocolos Modernos de Medición

## Preparación de datos

### Estructura de datos requerida

Los datos deben organizarse en una matriz sitios × especies donde:
- Filas = sitios/muestras
- Columnas = especies
- Valores = presencia/ausencia (0/1) o abundancias


``` r
# Cargar datos de ejemplo (usando datos incluidos en vegan)
data(dune)
data(dune.env)

# Visualizar estructura de los datos
print("Primeras 5 filas y 10 columnas de la matriz de especies:")
## [1] "Primeras 5 filas y 10 columnas de la matriz de especies:"
kable(dune[1:5, 1:10])
Achimill Agrostol Airaprae Alopgeni Anthodor Bellpere Bromhord Chenalbu Cirsarve Comapalu
1 0 0 0 0 0 0 0 0 0
3 0 0 2 0 3 4 0 0 0
0 4 0 7 0 2 0 0 0 0
0 8 0 2 0 2 3 0 2 0
2 0 0 0 4 2 2 0 0 0 
print("\nDimensiones de la matriz:")
## [1] "\nDimensiones de la matriz:"
print(paste("Sitios:", nrow(dune), "| Especies:", ncol(dune)))
## [1] "Sitios: 20 | Especies: 30"
print("\nInformación ambiental:")
## [1] "\nInformación ambiental:"
kable(head(dune.env))
A1 Moisture Management Use Manure
2.8 1 SF Haypastu 4
3.5 1 BF Haypastu 2
4.3 2 SF Haypastu 4
4.2 2 SF Haypastu 4
6.3 1 HF Hayfield 2
4.3 1 HF Haypastu 2

7.0.1 Control de calidad de datos

# Verificar datos faltantes
na_count <- sum(is.na(dune))
print(paste("Valores faltantes:", na_count))
## [1] "Valores faltantes: 0"
# Especies con baja frecuencia
species_freq <- colSums(dune > 0)
rare_species <- sum(species_freq <= 2)
print(paste("Especies raras (≤2 sitios):", rare_species))
## [1] "Especies raras (≤2 sitios): 5"
# Sitios vacíos
empty_sites <- sum(rowSums(dune) == 0)
print(paste("Sitios vacíos:", empty_sites))
## [1] "Sitios vacíos: 0"
# Visualizar distribución de riqueza por sitio
site_richness <- rowSums(dune > 0)

# Crear múltiples gráficos de diagnóstico
par(mfrow = c(2, 2))

# Histograma de riqueza
hist(site_richness, 
     main = "Distribución de riqueza por sitio", 
     xlab = "Número de especies", 
     ylab = "Frecuencia",
     col = "skyblue",
     breaks = 8)

# Distribución de frecuencias de especies
hist(species_freq,
     main = "Frecuencia de especies",
     xlab = "Número de sitios donde aparece",
     ylab = "Número de especies",
     col = "lightcoral",
     breaks = 8)

# Curva de acumulación de especies por sitio
plot(1:nrow(dune), cumsum(sort(site_richness, decreasing = TRUE)),
     type = "b", pch = 16,
     main = "Acumulación de riqueza",
     xlab = "Sitios (ordenados por riqueza)",
     ylab = "Riqueza acumulada",
     col = "darkgreen")

# Abundancia total vs riqueza
total_abundance <- rowSums(dune)
plot(site_richness, total_abundance,
     main = "Riqueza vs Abundancia",
     xlab = "Riqueza de especies",
     ylab = "Abundancia total",
     pch = 16, col = "purple")
abline(lm(total_abundance ~ site_richness), col = "red", lty = 2)

par(mfrow = c(1, 1))

7.1 Cálculo de índices básicos

# Calcular índices de disimilitud básicos
# Sørensen (equivalente a Bray-Curtis para datos binarios)
dune_pa <- decostand(dune, method = "pa")
sorensen_dist <- vegdist(dune_pa, method = "bray")
jaccard_dist <- vegdist(dune_pa, method = "jaccard")
bray_curtis_dist <- vegdist(dune, method = "bray")

# Crear resumen de distancias
dist_summary <- data.frame(
  Indice = c("Sørensen", "Jaccard", "Bray-Curtis"),
  Media = c(mean(sorensen_dist), mean(jaccard_dist), mean(bray_curtis_dist)),
  Min = c(min(sorensen_dist), min(jaccard_dist), min(bray_curtis_dist)),
  Max = c(max(sorensen_dist), max(jaccard_dist), max(bray_curtis_dist))
)

kable(dist_summary, digits = 3, 
      caption = "Resumen de índices de disimilitud")
Resumen de índices de disimilitud
Indice Media Min Max
Sørensen 0.552 0.083 1
Jaccard 0.688 0.154 1
Bray-Curtis 0.646 0.227 1

7.2 Interpretación de los Índices de Disimilitud

Los valores promedio de disimilitud se encuentran en el rango moderado a alto (0.552-0.688), indicando que las comunidades estudiadas presentan diferencias sustanciales en su composición de especies.

El índice de Sørensen muestra la menor disimilitud promedio (0.552), lo que refleja su mayor sensibilidad a las especies compartidas entre sitios. Este índice otorga mayor peso a las coincidencias, por lo que valores cercanos a 0.55 sugieren que, aunque existe diferenciación, las comunidades mantienen un núcleo de especies comunes. El valor mínimo de 0.083 indica que algunos pares de sitios son notablemente similares, mientras que el máximo de 1.0 revela la existencia de comunidades completamente diferentes sin especies compartidas.

El índice de Jaccard presenta la disimilitud más alta (0.688), resultado esperado dado que es más estricto al evaluar similitud, ya que no pondera doblemente las especies compartidas como lo hace Sørensen. La diferencia entre ambos índices (aproximadamente 0.14) es consistente con lo reportado en literatura ecológica y sugiere un nivel moderado de especies compartidas entre las comunidades.

Bray-Curtis, con un valor intermedio (0.646), incorpora información de abundancia además de presencia-ausencia, lo que proporciona una perspectiva más completa de la estructura comunitaria. Su valor mínimo más alto (0.227) comparado con los otros índices sugiere que, incluso en sitios con composiciones similares, existen diferencias notables en las abundancias relativas de las especies.

La amplitud de valores (rango mínimo-máximo) en los tres índices indica un gradiente continuo de diferenciación, desde comunidades muy similares hasta completamente diferentes. Esta variabilidad sugiere la presencia de factores ambientales o espaciales que generan heterogeneidad en el ensamblaje de especies, lo cual es fundamental para entender los procesos ecológicos subyacentes y tiene implicaciones importantes para estrategias de conservación y manejo del paisaje estudiado.

# Verificar que las variables de distancia existen
if(exists("sorensen_dist") && exists("jaccard_dist") && exists("bray_curtis_dist")) {
  
  # Crear data frame con las distancias
  dist_data <- data.frame(
    Sorensen = as.vector(sorensen_dist),
    Jaccard = as.vector(jaccard_dist),
    Bray_Curtis = as.vector(bray_curtis_dist)
  )
  
  # Transformar a formato largo para ggplot
  dist_long <- pivot_longer(dist_data, everything(), names_to = "Indice", values_to = "Distancia")
  
  # Visualizar distribuciones de distancias
  p1 <- ggplot(dist_long, aes(x = Distancia, fill = Indice)) +
    geom_histogram(alpha = 0.7, bins = 15, position = "identity") +
    facet_wrap(~Indice, scales = "free_y") +
    scale_fill_brewer(palette = "Set2") +
    theme_minimal() +
    labs(title = "Distribución de valores de disimilitud",
         x = "Valor de disimilitud", y = "Frecuencia") +
    theme(legend.position = "none")
  
  print(p1)
  
  # Matriz de correlaciones entre índices
  cor_matrix <- cor(dist_data)
  corrplot(cor_matrix, method = "color", type = "upper", 
           addCoef.col = "black", tl.cex = 0.8,
           title = "Correlaciones entre índices de disimilitud",
           mar = c(0,0,2,0))
  
} else {
  print("Error: Las variables de distancia no están disponibles.")
  print("Asegúrate de que el chunk 'basic-indices' se haya ejecutado correctamente.")
  
  # Mostrar qué variables existen
  available_vars <- c()
  if(exists("sorensen_dist")) available_vars <- c(available_vars, "sorensen_dist")
  if(exists("jaccard_dist")) available_vars <- c(available_vars, "jaccard_dist") 
  if(exists("bray_curtis_dist")) available_vars <- c(available_vars, "bray_curtis_dist")
  
  print(paste("Variables disponibles:", paste(available_vars, collapse = ", ")))
}

7.3 Análisis con el paquete betapart

# Convertir a presencia/ausencia
dune_pa <- decostand(dune, method = "pa")

# Calcular diversidad beta y sus componentes
beta_multi <- beta.multi(dune_pa, index.family = "sorensen")
beta_pair <- beta.pair(dune_pa, index.family = "sorensen")

print("Diversidad beta múltiple (Sørensen):")
## [1] "Diversidad beta múltiple (Sørensen):"
print(beta_multi)
## $beta.SIM
## [1] 0.8146504
## 
## $beta.SNE
## [1] 0.0421863
## 
## $beta.SOR
## [1] 0.8568367
# DIAGNÓSTICO: Ver qué contiene beta_multi
print("Estructura de beta_multi:")
## [1] "Estructura de beta_multi:"
print(str(beta_multi))
## List of 3
##  $ beta.SIM: num 0.815
##  $ beta.SNE: num 0.0422
##  $ beta.SOR: num 0.857
## NULL
# SOLUCIÓN: Extraer valores de forma segura
beta_total <- ifelse("beta.SOR" %in% names(beta_multi), beta_multi$beta.SOR, NA)
beta_sim <- ifelse("beta.SIM" %in% names(beta_multi), beta_multi$beta.SIM, NA)  
beta_nes <- ifelse("beta.NES" %in% names(beta_multi), beta_multi$beta.NES, NA)

# Crear data.frame solo con valores válidos
valores_validos <- c(beta_total, beta_sim, beta_nes)
componentes_nombres <- c("Beta total", "Recambio", "Anidamiento")

# Filtrar solo valores no NA
indices_validos <- !is.na(valores_validos)
beta_components <- data.frame(
  Componente = componentes_nombres[indices_validos],
  Valor = valores_validos[indices_validos]
)

# Verificar que tenemos datos para graficar
if(nrow(beta_components) > 0) {
  # Visualizar componentes de diversidad beta
  ggplot(beta_components, aes(x = Componente, y = Valor, fill = Componente)) +
    geom_col() +
    scale_fill_brewer(palette = "Set2") +
    theme_minimal() +
    labs(title = "Componentes de la Diversidad Beta",
         y = "Valor del índice") +
    theme(legend.position = "none")
} else {
  print("No hay componentes de diversidad beta para visualizar")
}

La diversidad beta total (~0.85) es muy alta, indicando diferencias sustanciales entre comunidades. Más revelador es que el componente de recambio (~0.85) representa prácticamente toda la diversidad beta, mientras que no hay componente de anidamiento visible. Esta dominancia del recambio sugiere que las diferencias entre sitios se deben principalmente a la sustitución de especies más que a diferencias en riqueza. Las comunidades tienen riquezas similares pero composiciones muy diferentes, indicando procesos ecológicos como especialización de hábitat, competencia interespecífica, o efectos de dispersión limitada. Este patrón es típico de ambientes con filtros ambientales fuertes que favorecen diferentes conjuntos de especies en cada sitio.

## Análisis de ordenación
# Análisis de Correspondencia (CA)
dune_ca <- cca(dune)

# Escalamiento multidimensional no métrico (NMDS)
dune_nmds <- metaMDS(dune, distance = "bray", k = 2)
## Run 0 stress 0.1192678 
## Run 1 stress 0.192224 
## Run 2 stress 0.1192678 
## ... Procrustes: rmse 5.026438e-05  max resid 0.0001571376 
## ... Similar to previous best
## Run 3 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 0.0202703  max resid 0.06496225 
## Run 4 stress 0.1192679 
## Run 5 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 6.697513e-06  max resid 2.000266e-05 
## ... Similar to previous best
## Run 6 stress 0.1192679 
## Run 7 stress 0.1192679 
## Run 8 stress 0.1192678 
## Run 9 stress 0.1183186 
## ... Procrustes: rmse 1.784713e-05  max resid 5.749849e-05 
## ... Similar to previous best
## Run 10 stress 0.1192679 
## Run 11 stress 0.1192678 
## Run 12 stress 0.3098501 
## Run 13 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 8.468079e-06  max resid 3.021793e-05 
## ... Similar to previous best
## Run 14 stress 0.1192678 
## Run 15 stress 0.1192678 
## Run 16 stress 0.1192678 
## Run 17 stress 0.1192678 
## Run 18 stress 0.2042114 
## Run 19 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 4.306659e-06  max resid 1.539611e-05 
## ... Similar to previous best
## Run 20 stress 0.1183186 
## ... Procrustes: rmse 2.933359e-06  max resid 1.039783e-05 
## ... Similar to previous best
## *** Best solution repeated 2 times
# Visualizar NMDS básico
plot(dune_nmds, type = "t", main = "NMDS - Sitios de dunas")

El análisis NMDS muestra la distribución de especies en el espacio multidimensional, revelando gradientes ecológicos subyacentes. La dispersión de especies como Airaprae, Empenigr, Hyporadi en la parte superior versus Chenalbu, Circarve en la inferior sugiere un gradiente ambiental que estructura las comunidades.

La distribución no aleatoria indica que factores ambientales específicos (posiblemente humedad, salinidad, o perturbación) determinan qué especies coexisten. Las especies agrupadas espacialmente probablemente comparten preferencias ecológicas similares.

# Crear gráfico mejorado con ggplot
nmds_scores <- scores(dune_nmds, display = "sites") %>%
  as.data.frame() %>%
  mutate(Site = rownames(.),
         Management = dune.env$Management)

ggplot(nmds_scores, aes(x = NMDS1, y = NMDS2, color = Management)) +
  geom_point(size = 3) +
  stat_ellipse(type = "norm", linetype = 2) +
  theme_minimal() +
  labs(title = "Ordenación NMDS de sitios de dunas",
       subtitle = paste("Stress =", round(dune_nmds$stress, 3))) +
  theme(legend.position = "bottom")
## Too few points to calculate an ellipse
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_path()`).

Este análisis es particularmente revelador sobre el efecto del manejo humano en la estructura comunitaria. El stress de 0.118 indica una representación excelente de las relaciones multidimensionales en 2D. Los patrones muestran:

  • NM (Natural Meadow) sitios agrupados principalmente en la derecha, sugiriendo comunidades distintivas en áreas no manejadas

  • SF (Standard Farming) y otros tipos de manejo forman grupos parcialmente superpuestos pero distinguibles

  • Las elipses de confianza muestran que diferentes intensidades de manejo crean nichos ecológicos específicos

  • La separación clara entre tipos de manejo confirma que las actividades humanas son un driver principal de la diferenciación comunitaria, más que factores puramente naturales.

7.4 Visualizaciones avanzadas

# 1. Heatmap de matriz de disimilitud
dist_matrix_full <- as.matrix(vegdist(dune, method = "bray"))

# Crear heatmap usando ggplot
dist_melted <- expand.grid(Site1 = rownames(dist_matrix_full), 
                          Site2 = colnames(dist_matrix_full)) %>%
  mutate(Distance = as.vector(dist_matrix_full))

ggplot(dist_melted, aes(x = Site1, y = Site2, fill = Distance)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "red", name = "Disimilitud\nBray-Curtis") +
  theme_minimal() +
  labs(title = "Matriz de disimilitud entre sitios",
       x = "Sitios", y = "Sitios") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text.y = element_text(angle = 0))

# 2. Red de similitud
library(igraph)
## 
## Attaching package: 'igraph'
## The following object is masked from 'package:tidyr':
## 
##     crossing
## The following objects are masked from 'package:dplyr':
## 
##     as_data_frame, groups, union
## The following object is masked from 'package:vegan':
## 
##     diversity
## The following object is masked from 'package:permute':
## 
##     permute
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
# Convertir distancias a similitudes
similarity_matrix <- 1 - dist_matrix_full
# Mantener solo conexiones fuertes (similitud > 0.6)
similarity_matrix[similarity_matrix < 0.6] <- 0
diag(similarity_matrix) <- 0

# Crear grafo
graph_net <- graph_from_adjacency_matrix(similarity_matrix, 
                                        weighted = TRUE, 
                                        mode = "undirected")

# Configurar colores por tipo de manejo
vertex_colors <- rainbow(length(unique(dune.env$Management)))[as.numeric(factor(dune.env$Management))]

plot(graph_net,
     vertex.color = vertex_colors,
     vertex.size = 8,
     vertex.label = 1:nrow(dune),
     edge.width = E(graph_net)$weight * 3,
     main = "Red de similitud entre sitios (>60% similitud)",
     layout = layout_with_fr)

legend("topright", 
       legend = levels(factor(dune.env$Management)),
       col = rainbow(length(unique(dune.env$Management))),
       pch = 19, cex = 0.8)

La red revela patrones de conectividad ecológica fascinantes. Solo se muestran conexiones con >60% similitud, indicando que la mayoría de sitios son bastante diferentes entre sí. Observaciones clave:

  • Existen clusters de sitios similares conectados densamente
  • Algunos sitios actúan como “puentes” conectando diferentes grupos
  • Sitios aislados (sin conexiones) representan comunidades únicas
  • La distribución por colores de manejo muestra que sitios con el mismo tipo de manejo tienden a conectarse más

Esta estructura sugiere metacomunidades donde ciertos sitios funcionan como fuentes de colonización para otros, y donde el tipo de manejo determina la compatibilidad ecológica entre sitios.

# 3. Análisis de agrupamiento con bootstrapping
# SOLUCIÓN: Usar método de distancia válido para pvclust
library(pvclust)

# Opción 1: Usar distancia euclidiana (válida para pvclust)
dune_clust_euclidean <- pvclust(t(dune), method.hclust = "average", 
                               method.dist = "euclidean", nboot = 100)
## Bootstrap (r = 0.5)... Done.
## Bootstrap (r = 0.6)... Done.
## Bootstrap (r = 0.7)... Done.
## Bootstrap (r = 0.8)... Done.
## Bootstrap (r = 0.9)... Done.
## Bootstrap (r = 1.0)... Done.
## Bootstrap (r = 1.1)... Done.
## Bootstrap (r = 1.2)... Done.
## Bootstrap (r = 1.3)... Done.
## Bootstrap (r = 1.4)... Done.
plot(dune_clust_euclidean, main = "Dendrograma con soporte bootstrap (Euclidiana)")
pvrect(dune_clust_euclidean, alpha = 0.95)

El dendrograma muestra relaciones jerárquicas entre sitios con soporte estadístico robusto. Los valores de bootstrap (en rojo y verde) indican la confianza en cada agrupamiento:

  • Valores altos (>80) en varios nodos confirman agrupamientos estadísticamente sólidos
  • Los rectángulos grises delimitan grupos principales con alta confianza
  • La altura variable de las ramificaciones indica diferentes grados de disimilitud

El uso de distancia euclidiana en lugar de índices ecológicos puede explicar algunos agrupamientos inesperados, pero la estructura general sugiere procesos históricos o gradientes ambientales que han moldeado patrones jerárquicos de similitud entre comunidades.

# Opción 2: Clustering manual con Bray-Curtis y bootstrap
# Crear función para bootstrap manual
bootstrap_cluster <- function(data, nboot = 100) {
  n_sites <- nrow(data)
  dist_bray <- vegdist(data, method = "bray")
  hclust_result <- hclust(dist_bray, method = "average")
  
  # Bootstrap simple
  bootstrap_supports <- numeric(length(hclust_result$height))
  
  for(i in 1:nboot) {
    # Resample especies
    boot_indices <- sample(ncol(data), ncol(data), replace = TRUE)
    boot_data <- data[, boot_indices]
    
    # Calcular clustering bootstrap
    boot_dist <- vegdist(boot_data, method = "bray")
    boot_hclust <- hclust(boot_dist, method = "average")
    
    # Comparar topologías (simplificado)
    # En un análisis real, se compararían todas las particiones
  }
  
  return(hclust_result)
}

# Ejecutar clustering manual
hclust_bray <- hclust(vegdist(dune, method = "bray"), method = "average")
plot(hclust_bray, main = "Dendrograma (Bray-Curtis, sin bootstrap)",
     xlab = "Sitios", ylab = "Distancia")

# Agregar rectángulos para grupos
rect.hclust(hclust_bray, k = 4, border = "red")

8 Ejemplos Prácticos Avanzados

8.1 Ejemplo 1: Diversidad beta a lo largo de gradientes ambientales

# Crear gradiente ambiental sintético
set.seed(42)
n_sites <- 20
n_species <- 15

# Simular gradiente de pH
ph_gradient <- seq(4, 8, length.out = n_sites)

# Simular respuestas de especies al pH
species_optima <- runif(n_species, min = 4.5, max = 7.5)
species_tolerance <- runif(n_species, min = 0.5, max = 1.5)

# Generar matriz de comunidades
comm_gradient <- matrix(0, nrow = n_sites, ncol = n_species)

for(i in 1:n_sites) {
  for(j in 1:n_species) {
    prob <- exp(-((ph_gradient[i] - species_optima[j])^2) / (2 * species_tolerance[j]^2))
    comm_gradient[i, j] <- rbinom(1, 1, prob * 0.8)  # probabilidad máxima de 0.8
  }
}

rownames(comm_gradient) <- paste("Site", 1:n_sites)
colnames(comm_gradient) <- paste("Sp", 1:n_species)

# Calcular diversidad beta a lo largo del gradiente
beta_gradient <- beta.multi(comm_gradient)

# Calcular distancias por pares y correlacionar con diferencias de pH
dist_beta <- as.matrix(beta.pair(comm_gradient, index.family = "sorensen")$beta.sor)
ph_diff <- as.matrix(dist(ph_gradient))

# Correlación entre distancia ambiental y beta diversidad
correlation <- cor(as.vector(ph_diff[upper.tri(ph_diff)]), 
                  as.vector(dist_beta[upper.tri(dist_beta)]))

print(paste("Correlación pH vs Beta diversidad:", round(correlation, 3)))
## [1] "Correlación pH vs Beta diversidad: 0.562"
# Visualizar
plot_data <- data.frame(
  pH_diff = as.vector(ph_diff[upper.tri(ph_diff)]),
  Beta_div = as.vector(dist_beta[upper.tri(dist_beta)])
)

ggplot(plot_data, aes(x = pH_diff, y = Beta_div)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_minimal() +
  labs(x = "Diferencia de pH entre sitios",
       y = "Diversidad Beta (Sørensen)",
       title = "Relación entre gradiente ambiental y diversidad beta") +
  annotate("text", x = max(plot_data$pH_diff) * 0.7, y = max(plot_data$Beta_div) * 0.9,
           label = paste("r =", round(correlation, 3)))
## `geom_smooth()` using formula = 'y ~ x'

8.2 Ejemplo 2: Partición de diversidad beta en escalas espaciales

# Simular datos espaciales jerárquicos
set.seed(100)
n_regions <- 3
sites_per_region <- 6
n_species <- 20

# Crear estructura jerárquica
region_labels <- rep(1:n_regions, each = sites_per_region)
site_labels <- paste("R", rep(1:n_regions, each = sites_per_region), 
                    "S", rep(1:sites_per_region, n_regions), sep = "")

# Pool de especies por región
regional_pools <- list(
  region1 = 1:12,
  region2 = 8:20,
  region3 = c(1:5, 15:20)
)

# Generar comunidades
comm_hierarchical <- matrix(0, nrow = length(site_labels), ncol = n_species)
rownames(comm_hierarchical) <- site_labels
colnames(comm_hierarchical) <- paste("Sp", 1:n_species)

for(i in 1:length(site_labels)) {
  region <- region_labels[i]
  available_species <- regional_pools[[region]]
  n_local <- sample(4:8, 1)  # riqueza local variable
  local_species <- sample(available_species, n_local, replace = FALSE)
  comm_hierarchical[i, local_species] <- 1
}

# Calcular partición aditiva de diversidad
alpha_mean <- mean(rowSums(comm_hierarchical))
gamma_total <- sum(colSums(comm_hierarchical) > 0)
beta_total <- gamma_total - alpha_mean

# Por regiones
beta_within_regions <- numeric(n_regions)
for(r in 1:n_regions) {
  sites_in_region <- which(region_labels == r)
  comm_region <- comm_hierarchical[sites_in_region, ]
  alpha_region <- mean(rowSums(comm_region))
  gamma_region <- sum(colSums(comm_region) > 0)
  beta_within_regions[r] <- gamma_region - alpha_region
}

beta_within_mean <- mean(beta_within_regions)
beta_among_regions <- beta_total - beta_within_mean

# Crear resumen
partition_summary <- data.frame(
  Componente = c("Alpha promedio", "Beta dentro regiones", "Beta entre regiones", "Gamma total"),
  Valor = c(alpha_mean, beta_within_mean, beta_among_regions, gamma_total),
  Proporcion = c(alpha_mean/gamma_total, beta_within_mean/gamma_total, 
                beta_among_regions/gamma_total, 1)
)

kable(partition_summary, digits = 2, 
      caption = "Partición aditiva de la diversidad")
Partición aditiva de la diversidad
Componente Valor Proporcion
Alpha promedio 5.61 0.28
Beta dentro regiones 6.39 0.32
Beta entre regiones 8.00 0.40
Gamma total 20.00 1.00 
# Visualizar
ggplot(partition_summary[-4,], aes(x = Componente, y = Proporcion, fill = Componente)) +
  geom_col() +
  scale_fill_brewer(palette = "Pastel1") +
  theme_minimal() +
  labs(title = "Partición de la diversidad en escalas espaciales",
       y = "Proporción de la diversidad gamma") +
  theme(legend.position = "none", 
        axis.text.x = element_text(angle = 45, hjust = 1))

8.3 Ejemplo 3: Efectos de la distancia y el ambiente

# Datos del paquete vegan: mites (ácaros del suelo)
data(mite)
data(mite.env)
data(mite.xy)

# Preparar datos
mite_pa <- decostand(mite, method = "pa")

# Calcular matrices de distancia
dist_species <- vegdist(mite_pa, method = "jaccard")
dist_env <- vegdist(scale(mite.env[,1:2]), method = "euclidean")  # usar solo variables numéricas
dist_geo <- dist(mite.xy)

# Prueba de Mantel
mantel_env <- mantel(dist_species, dist_env, method = "pearson", permutations = 999)
mantel_geo <- mantel(dist_species, dist_geo, method = "pearson", permutations = 999)

print("Prueba de Mantel - Ambiente:")
## [1] "Prueba de Mantel - Ambiente:"
print(mantel_env)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = dist_species, ydis = dist_env, method = "pearson",      permutations = 999) 
## 
## Mantel statistic r: 0.4263 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0597 0.0773 0.1005 0.1200 
## Permutation: free
## Number of permutations: 999
print("\nPrueba de Mantel - Geografía:")
## [1] "\nPrueba de Mantel - Geografía:"
print(mantel_geo)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = dist_species, ydis = dist_geo, method = "pearson",      permutations = 999) 
## 
## Mantel statistic r: 0.6287 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0443 0.0604 0.0754 0.0961 
## Permutation: free
## Number of permutations: 999
# Mantel parcial (ambiente controlando por geografía)
mantel_partial <- mantel.partial(dist_species, dist_env, dist_geo, 
                                method = "pearson", permutations = 999)

print("\nMantel parcial (Ambiente | Geografía):")
## [1] "\nMantel parcial (Ambiente | Geografía):"
print(mantel_partial)
## 
## Partial Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel.partial(xdis = dist_species, ydis = dist_env, zdis = dist_geo,      method = "pearson", permutations = 999) 
## 
## Mantel statistic r: 0.3323 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0624 0.0793 0.0975 0.1121 
## Permutation: free
## Number of permutations: 999
# Visualizar relaciones
par(mfrow = c(1, 2))
plot(dist_env, dist_species, pch = 16, col = rgb(0, 0, 1, 0.6),
     xlab = "Distancia ambiental", ylab = "Distancia en composición de especies",
     main = paste("Ambiente (r =", round(mantel_env$statistic, 3), ")"))
abline(lm(as.vector(dist_species) ~ as.vector(dist_env)), col = "red")

plot(dist_geo, dist_species, pch = 16, col = rgb(1, 0, 0, 0.6),
     xlab = "Distancia geográfica", ylab = "Distancia en composición de especies",
     main = paste("Geografía (r =", round(mantel_geo$statistic, 3), ")"))
abline(lm(as.vector(dist_species) ~ as.vector(dist_geo)), col = "blue")

9 Ejercicios Prácticos

9.1 Ejercicio 1: Análisis básico de diversidad beta

Objetivo: Calcular y comparar diferentes índices de diversidad beta.

# INSTRUCCIONES:
# 1. Carga el conjunto de datos 'varespec' del paquete vegan
# 2. Calcula las matrices de distancia usando Jaccard, Bray-Curtis y Sørensen
# 3. Compara las correlaciones entre estos índices
# 4. Crea un dendrograma usando una de las matrices de distancia

# Tu código aquí:
data(varespec)

# Paso 1: Explorar los datos
# ¿Cuántos sitios y especies hay?
# ¿Hay sitios vacíos o especies muy raras?

# Paso 2: Calcular índices de disimilitud
# jaccard_dist <- ?
# bray_dist <- ?
# sorensen_dist <- ?

# Paso 3: Matriz de correlaciones
# dist_matrix <- cbind(...)
# cor_matrix <- cor(dist_matrix)
# corrplot(cor_matrix)

# Paso 4: Dendrograma
# hclust_result <- ?
# plot(hclust_result)

9.1.1 Solución del Ejercicio 1

# Cargar datos
data(varespec)

# Explorar datos
print(paste("Sitios:", nrow(varespec), "| Especies:", ncol(varespec)))
## [1] "Sitios: 24 | Especies: 44"
print(paste("Sitios vacíos:", sum(rowSums(varespec) == 0)))
## [1] "Sitios vacíos: 0"
print(paste("Especies en un solo sitio:", sum(colSums(varespec > 0) == 1)))
## [1] "Especies en un solo sitio: 0"
# Calcular índices
varespec_pa <- decostand(varespec, method = "pa")
jaccard_dist <- vegdist(varespec_pa, method = "jaccard")
bray_dist <- vegdist(varespec, method = "bray")

# Para Sørensen con datos de abundancia, transformamos primero
sorensen_dist <- vegdist(varespec_pa, method = "bray")  # equivalente para datos binarios

# Matriz de correlaciones
dist_matrix <- cbind(
  Jaccard = as.vector(jaccard_dist),
  Bray_Curtis = as.vector(bray_dist),
  Sorensen = as.vector(sorensen_dist)
)

cor_matrix <- cor(dist_matrix)
print("Matriz de correlaciones:")
## [1] "Matriz de correlaciones:"
kable(cor_matrix, digits = 3)
Jaccard Bray_Curtis Sorensen
Jaccard 1.000 0.221 0.997
Bray_Curtis 0.221 1.000 0.220
Sorensen 0.997 0.220 1.000 
# Visualizar correlaciones
corrplot(cor_matrix, method = "color", type = "upper", 
         addCoef.col = "black", tl.cex = 0.8)

# Dendrograma
hclust_result <- hclust(bray_dist, method = "average")
plot(hclust_result, main = "Dendrograma (Bray-Curtis, UPGMA)",
     xlab = "Sitios", ylab = "Distancia")

9.2 Ejercicio 2: Descomposición de diversidad beta

Objetivo: Usar el paquete betapart para descomponer la diversidad beta.

# INSTRUCCIONES:
# 1. Usa los datos 'dune' para calcular la diversidad beta múltiple
# 2. Descompón la diversidad beta en sus componentes de recambio y anidamiento
# 3. Identifica qué pares de sitios tienen mayor recambio vs. anidamiento
# 4. Relaciona estos patrones con las variables ambientales

# Tu código aquí:
data(dune)
data(dune.env)

# Paso 1: Diversidad beta múltiple
# beta_result <- ?

# Paso 2: Diversidad beta por pares
# beta_pairs <- ?

# Paso 3: Identificar patrones extremos
# ¿Qué pares tienen mayor recambio?
# ¿Qué pares tienen mayor anidamiento?

# Paso 4: Relación con ambiente
# Crear matriz de distancias ambientales y compararla

9.2.1 Solución del Ejercicio 2

# Convertir a presencia/ausencia
dune_pa <- decostand(dune, method = "pa")

# Diversidad beta múltiple
beta_multi <- beta.multi(dune_pa, index.family = "sorensen")
print("Diversidad beta múltiple:")
## [1] "Diversidad beta múltiple:"
print(beta_multi)
## $beta.SIM
## [1] 0.8146504
## 
## $beta.SNE
## [1] 0.0421863
## 
## $beta.SOR
## [1] 0.8568367
# Diversidad beta por pares
beta_pairs <- beta.pair(dune_pa, index.family = "sorensen")

# DIAGNÓSTICO: Verificar qué componentes están disponibles
print("Componentes disponibles en beta_pairs:")
## [1] "Componentes disponibles en beta_pairs:"
print(names(beta_pairs))
## [1] "beta.sim" "beta.sne" "beta.sor"
print("Estructura de beta_pairs:")
## [1] "Estructura de beta_pairs:"
print(str(beta_pairs))
## List of 3
##  $ beta.sim: 'dist' Named num [1:190] 0 0.2 0.2 0 0.2 0.2 0.4 0.2 0.2 0.6 ...
##   ..- attr(*, "Labels")= chr [1:20] "1" "2" "3" "4" ...
##   ..- attr(*, "Size")= int 20
##   ..- attr(*, "call")= language as.dist.default(m = beta.sim)
##   ..- attr(*, "Diag")= logi FALSE
##   ..- attr(*, "Upper")= logi FALSE
##  $ beta.sne: 'dist' Named num [1:190] 0.333 0.267 0.356 0.474 0.3 ...
##   ..- attr(*, "Labels")= chr [1:20] "1" "2" "3" "4" ...
##   ..- attr(*, "Size")= int 20
##   ..- attr(*, "call")= language as.dist.default(m = beta.sne)
##   ..- attr(*, "Diag")= logi FALSE
##   ..- attr(*, "Upper")= logi FALSE
##  $ beta.sor: 'dist' Named num [1:190] 0.333 0.467 0.556 0.474 0.5 ...
##   ..- attr(*, "Labels")= chr [1:20] "1" "2" "3" "4" ...
##   ..- attr(*, "Size")= int 20
##   ..- attr(*, "call")= language as.dist.default(m = beta.sor)
##   ..- attr(*, "Diag")= logi FALSE
##   ..- attr(*, "Upper")= logi FALSE
## NULL
# SOLUCIÓN: Verificar componentes antes de usarlos
if("beta.sor" %in% names(beta_pairs) && !is.null(beta_pairs$beta.sor)) {
  sor_matrix <- as.matrix(beta_pairs$beta.sor)
} else {
  print("Componente beta.sor no disponible")
  sor_matrix <- NULL
}

if("beta.sim" %in% names(beta_pairs) && !is.null(beta_pairs$beta.sim)) {
  sim_matrix <- as.matrix(beta_pairs$beta.sim)
} else {
  print("Componente beta.sim no disponible")
  sim_matrix <- NULL
}

if("beta.nes" %in% names(beta_pairs) && !is.null(beta_pairs$beta.nes)) {
  nes_matrix <- as.matrix(beta_pairs$beta.nes)
} else {
  print("Componente beta.nes no disponible")
  nes_matrix <- NULL
}
## [1] "Componente beta.nes no disponible"
# Continuar solo si tenemos las matrices necesarias
if(!is.null(sor_matrix) && !is.null(sim_matrix)) {
  # Identificar pares extremos
  # Encontrar máximos (excluyendo diagonal)
  diag(sor_matrix) <- NA
  diag(sim_matrix) <- NA
  if(!is.null(nes_matrix)) {
    diag(nes_matrix) <- NA
  }

  max_turnover <- which(sim_matrix == max(sim_matrix, na.rm = TRUE), arr.ind = TRUE)[1,]
  
  print(paste("Máximo recambio entre sitios:", 
             rownames(dune)[max_turnover[1]], "y", rownames(dune)[max_turnover[2]]))
  print(paste("Valor:", round(max(sim_matrix, na.rm = TRUE), 3)))

  # Solo mostrar anidamiento si está disponible
  if(!is.null(nes_matrix)) {
    max_nestedness <- which(nes_matrix == max(nes_matrix, na.rm = TRUE), arr.ind = TRUE)[1,]
    print(paste("Máximo anidamiento entre sitios:", 
               rownames(dune)[max_nestedness[1]], "y", rownames(dune)[max_nestedness[2]]))
    print(paste("Valor:", round(max(nes_matrix, na.rm = TRUE), 3)))
  } else {
    print("Componente de anidamiento no disponible para este análisis")
  }

  # Relación con variables ambientales
  # Crear matriz de distancia ambiental usando variables categóricas
  dune_env_numeric <- model.matrix(~ Management + Use + Manure - 1, data = dune.env)
  env_dist <- vegdist(dune_env_numeric, method = "euclidean")

  # DIAGNÓSTICO: Verificar dimensiones
  print("Diagnóstico de dimensiones:")
  print(paste("Dimensiones sor_matrix:", nrow(sor_matrix), "x", ncol(sor_matrix)))
  print(paste("Dimensiones sim_matrix:", nrow(sim_matrix), "x", ncol(sim_matrix)))
  print(paste("Longitud env_dist:", length(env_dist)))
  print(paste("Sitios en dune:", nrow(dune)))
  print(paste("Sitios en dune.env:", nrow(dune.env)))

  # SOLUCIÓN: Usar solo triángulo superior para correlaciones
  # Convertir matrices de distancia a formato vectorial consistente
  if(nrow(sor_matrix) == nrow(dune) && length(env_dist) == choose(nrow(dune), 2)) {
    # Usar triángulo superior de las matrices (excluyendo diagonal)
    sor_vector <- sor_matrix[upper.tri(sor_matrix)]
    sim_vector <- sim_matrix[upper.tri(sim_matrix)]
    env_vector <- as.vector(env_dist)
    
    # Verificar que tengan la misma longitud
    if(length(sor_vector) == length(env_vector)) {
      # Correlación con componentes de beta diversidad
      cor_total <- cor(sor_vector, env_vector, use = "complete.obs")
      cor_turnover <- cor(sim_vector, env_vector, use = "complete.obs")
      
      if(!is.null(nes_matrix)) {
        nes_vector <- nes_matrix[upper.tri(nes_matrix)]
        cor_nestedness <- cor(nes_vector, env_vector, use = "complete.obs")
        correlations <- data.frame(
          Componente = c("Beta total", "Recambio", "Anidamiento"),
          Correlacion_ambiente = c(cor_total, cor_turnover, cor_nestedness)
        )
      } else {
        correlations <- data.frame(
          Componente = c("Beta total", "Recambio"),
          Correlacion_ambiente = c(cor_total, cor_turnover)
        )
      }
    } else {
      print("Error: Las longitudes de los vectores no coinciden")
      print(paste("Longitud beta vectors:", length(sor_vector)))
      print(paste("Longitud env vector:", length(env_vector)))
      
      # Análisis alternativo simplificado
      correlations <- data.frame(
        Componente = c("Error de dimensiones"),
        Correlacion_ambiente = c(NA)
      )
    }
  } else {
    print("Error: Dimensiones incompatibles entre matrices")
    # Crear matrices de distancia directamente desde los datos originales
    beta_dist_direct <- vegdist(dune_pa, method = "jaccard")
    env_dist_direct <- vegdist(dune_env_numeric, method = "euclidean")
    
    if(length(beta_dist_direct) == length(env_dist_direct)) {
      cor_direct <- cor(as.vector(beta_dist_direct), as.vector(env_dist_direct), use = "complete.obs")
      correlations <- data.frame(
        Componente = c("Beta diversidad (directo)"),
        Correlacion_ambiente = c(cor_direct)
      )
    } else {
      correlations <- data.frame(
        Componente = c("Error dimensional"),
        Correlacion_ambiente = c(NA)
      )
    }
  }

  kable(correlations, digits = 3, 
        caption = "Correlación entre componentes de beta diversidad y ambiente")

} else {
  print("Error: No se pudieron obtener las matrices de diversidad beta necesarias")
  
  # Análisis alternativo usando funciones básicas
  print("Realizando análisis alternativo...")
  
  # Calcular diversidad beta usando vegdist directamente
  beta_basic <- vegdist(dune_pa, method = "jaccard")
  
  print("Resumen de diversidad beta (Jaccard):")
  print(summary(beta_basic))
  
  # Relación con ambiente usando análisis básico
  dune_env_numeric <- model.matrix(~ Management + Use + Manure - 1, data = dune.env)
  env_dist <- vegdist(dune_env_numeric, method = "euclidean")
  
  cor_basic <- cor(as.vector(beta_basic), as.vector(env_dist))
  print(paste("Correlación básica beta-ambiente:", round(cor_basic, 3)))
}
## [1] "Máximo recambio entre sitios: 14 y 1"
## [1] "Valor: 1"
## [1] "Componente de anidamiento no disponible para este análisis"
## [1] "Diagnóstico de dimensiones:"
## [1] "Dimensiones sor_matrix: 20 x 20"
## [1] "Dimensiones sim_matrix: 20 x 20"
## [1] "Longitud env_dist: 190"
## [1] "Sitios en dune: 20"
## [1] "Sitios en dune.env: 20"
Correlación entre componentes de beta diversidad y ambiente
Componente Correlacion_ambiente
Beta total 0.047
Recambio 0.033

9.3 Ejercicio 3: Análisis de la diversidad beta temporal

Objetivo: Simular y analizar cambios temporales en diversidad beta.

# INSTRUCCIONES:
# 1. Simula una serie temporal de comunidades (ej. 10 años, 5 sitios, 15 especies)
# 2. Incluye diferentes tipos de cambio: recambio gradual, perturbaciones, recuperación
# 3. Calcula la diversidad beta entre años consecutivos
# 4. Visualiza los patrones temporales

# Tu código aquí:
set.seed(123)
n_years <- 10
n_sites <- 5
n_species <- 15

# Paso 1: Crear serie temporal
# ¿Cómo simular diferentes escenarios de cambio?

# Paso 2: Calcular beta diversidad temporal
# Para cada par de años consecutivos

# Paso 3: Visualizar patrones
# Gráfico de líneas mostrando cambio temporal en beta diversidad

9.3.1 Solución del Ejercicio 3

set.seed(123)
n_years <- 10
n_sites <- 5
n_species <- 15

# Simular serie temporal con diferentes tipos de cambio
temporal_communities <- array(0, dim = c(n_years, n_sites, n_species))

# Año inicial
for(i in 1:n_sites) {
  n_species_initial <- sample(5:10, 1)
  species_present <- sample(1:n_species, n_species_initial)
  temporal_communities[1, i, species_present] <- 1
}

# Evolución temporal
for(year in 2:n_years) {
  for(site in 1:n_sites) {
    prev_community <- temporal_communities[year-1, site, ]
    
    # Diferentes tipos de cambio según el año
    if(year <= 4) {
      # Cambio gradual (recambio lento)
      prob_extinction <- 0.1
      prob_colonization <- 0.1
    } else if(year == 5) {
      # Perturbación (mayor extinción)
      prob_extinction <- 0.4
      prob_colonization <- 0.1
    } else {
      # Recuperación (mayor colonización)
      prob_extinction <- 0.1
      prob_colonization <- 0.3
    }
    
    # Aplicar cambios
    new_community <- prev_community
    
    # Extinción
    present_species <- which(prev_community == 1)
    extinctions <- rbinom(length(present_species), 1, prob_extinction)
    new_community[present_species[extinctions == 1]] <- 0
    
    # Colonización
    absent_species <- which(prev_community == 0)
    colonizations <- rbinom(length(absent_species), 1, prob_colonization)
    new_community[absent_species[colonizations == 1]] <- 1
    
    temporal_communities[year, site, ] <- new_community
  }
}

# Calcular beta diversidad temporal
beta_temporal <- numeric(n_years - 1)
turnover_temporal <- numeric(n_years - 1)
nestedness_temporal <- numeric(n_years - 1)

for(year in 2:n_years) {
  # Combinar sitios de años consecutivos
  comm_t1 <- temporal_communities[year-1, , ]
  comm_t2 <- temporal_communities[year, , ]
  
  # Combinar en una matriz para beta.multi
  combined_comm <- rbind(as.data.frame(comm_t1), as.data.frame(comm_t2))
  
  # SOLUCIÓN: Verificar que hay datos válidos antes de calcular
  if(nrow(combined_comm) > 1 && ncol(combined_comm) > 0 && sum(combined_comm) > 0) {
    # Intentar calcular beta diversidad
    tryCatch({
      beta_result <- beta.multi(combined_comm, index.family = "sorensen")
      
      # DIAGNÓSTICO: Verificar qué componentes están disponibles
      if(year == 2) {  # Solo imprimir una vez para no saturar la salida
        print("Componentes disponibles en beta_result:")
        print(names(beta_result))
      }
      
      # Asignar valores de forma segura
      beta_temporal[year-1] <- if("beta.SOR" %in% names(beta_result) && !is.null(beta_result$beta.SOR)) {
        beta_result$beta.SOR
      } else {
        NA
      }
      
      turnover_temporal[year-1] <- if("beta.SIM" %in% names(beta_result) && !is.null(beta_result$beta.SIM)) {
        beta_result$beta.SIM
      } else {
        NA
      }
      
      nestedness_temporal[year-1] <- if("beta.NES" %in% names(beta_result) && !is.null(beta_result$beta.NES) && length(beta_result$beta.NES) > 0) {
        beta_result$beta.NES
      } else {
        NA
      }
      
    }, error = function(e) {
      # Si beta.multi falla, calcular índices básicos
      if(year == 2) {
        print(paste("Error en beta.multi para año", year, ":", e$message))
        print("Usando métodos alternativos...")
      }
      
      # Calcular diversidad beta básica usando vegdist
      basic_beta <- vegdist(combined_comm, method = "jaccard")
      beta_temporal[year-1] <<- ifelse(length(basic_beta) > 0, mean(basic_beta), NA)
      turnover_temporal[year-1] <<- NA  # No se puede calcular sin betapart
      nestedness_temporal[year-1] <<- NA  # No se puede calcular sin betapart
    })
  } else {
    # No hay datos válidos para calcular beta diversidad
    beta_temporal[year-1] <- NA
    turnover_temporal[year-1] <- NA
    nestedness_temporal[year-1] <- NA
  }
}
## [1] "Componentes disponibles en beta_result:"
## [1] "beta.SIM" "beta.SNE" "beta.SOR"
# Crear dataframe para visualización (solo con valores no NA)
temporal_data <- data.frame(
  Year_transition = paste(1:(n_years-1), "to", 2:n_years),
  Year_numeric = 1:(n_years-1),
  Beta_total = beta_temporal,
  Turnover = turnover_temporal,
  Nestedness = nestedness_temporal
)

# Filtrar filas con datos válidos
temporal_data_clean <- temporal_data[!is.na(temporal_data$Beta_total), ]

if(nrow(temporal_data_clean) > 0) {
  # Transformar a formato largo para ggplot (solo componentes disponibles)
  cols_to_plot <- c("Beta_total")
  if(any(!is.na(temporal_data_clean$Turnover))) {
    cols_to_plot <- c(cols_to_plot, "Turnover")
  }
  if(any(!is.na(temporal_data_clean$Nestedness))) {
    cols_to_plot <- c(cols_to_plot, "Nestedness")
  }
  
  temporal_long <- temporal_data_clean %>%
    select(Year_numeric, all_of(cols_to_plot)) %>%
    pivot_longer(cols = -Year_numeric, names_to = "Component", values_to = "Value") %>%
    filter(!is.na(Value))

  if(nrow(temporal_long) > 0) {
    # Visualizar
    ggplot(temporal_long, aes(x = Year_numeric, y = Value, color = Component)) +
      geom_line(size = 1.2) +
      geom_point(size = 3) +
      scale_color_brewer(palette = "Dark2") +
      theme_minimal() +
      labs(x = "Transición temporal", y = "Diversidad Beta",
           title = "Diversidad Beta Temporal",
           subtitle = "Cambios en composición de especies a través del tiempo") +
      annotate("rect", xmin = 3.5, xmax = 4.5, ymin = -Inf, ymax = Inf,
               alpha = 0.2, fill = "red") +
      annotate("text", x = 4, y = max(temporal_long$Value, na.rm = TRUE) * 0.9, 
               label = "Perturbación", hjust = 0.5)
  } else {
    print("No hay datos válidos para visualizar")
    plot.new()
    text(0.5, 0.5, "No hay datos válidos para la visualización", cex = 1.5)
  }
} else {
  print("No se pudieron calcular valores de diversidad beta temporal")
  plot.new()
  text(0.5, 0.5, "Error: No se pudieron calcular valores de diversidad beta", cex = 1.5)
}
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Tabla resumen (solo con datos válidos)
print("Resumen de diversidad beta temporal:")
## [1] "Resumen de diversidad beta temporal:"
kable(temporal_data, digits = 3, 
      caption = "Diversidad beta temporal por componentes (NA = no calculable)")
Diversidad beta temporal por componentes (NA = no calculable)
Year_transition Year_numeric Beta_total Turnover Nestedness
1 to 2 1 0.692 0.618 NA
2 to 3 2 0.692 0.648 NA
3 to 4 3 0.671 0.627 NA
4 to 5 4 0.705 0.601 NA
5 to 6 5 0.734 0.655 NA
6 to 7 6 0.683 0.619 NA
7 to 8 7 0.618 0.531 NA
8 to 9 8 0.533 0.405 NA
9 to 10 9 0.466 0.280 NA 
# Estadísticas básicas
valid_beta <- temporal_data$Beta_total[!is.na(temporal_data$Beta_total)]
if(length(valid_beta) > 0) {
  print(paste("Diversidad beta promedio:", round(mean(valid_beta), 3)))
  print(paste("Máxima diversidad beta:", round(max(valid_beta), 3)))
  print(paste("Número de transiciones válidas:", length(valid_beta)))
} else {
  print("No se calcularon valores válidos de diversidad beta")
}
## [1] "Diversidad beta promedio: 0.644"
## [1] "Máxima diversidad beta: 0.734"
## [1] "Número de transiciones válidas: 9"

10 Interpretación y Aplicaciones

10.1 Criterios de interpretación

10.1.1 Valores de diversidad beta

  • β < 0.2: Baja diferenciación (comunidades similares)
  • 0.2 ≤ β < 0.5: Diferenciación moderada
  • 0.5 ≤ β < 0.8: Alta diferenciación
  • β ≥ 0.8: Diferenciación muy alta (comunidades muy distintas)

10.1.2 Dominancia de componentes

  • Recambio dominante: Sugiere procesos de sustitución de especies
  • Anidamiento dominante: Indica diferencias en riqueza más que sustitución

10.2 Aplicaciones en conservación

# Ejemplo: Priorización de sitios para conservación
set.seed(456)

# Simular comunidades en diferentes fragmentos
n_fragments <- 12
n_species_pool <- 25

# Crear comunidades con diferentes niveles de fragmentación
communities_conservation <- matrix(0, nrow = n_fragments, ncol = n_species_pool)

# Fragmentos grandes (sitios 1-4): alta diversidad, especies especialistas
for(i in 1:4) {
  n_local <- sample(12:18, 1)
  specialists <- sample(1:10, min(n_local, 8))  # especies especialistas
  generalists <- sample(11:n_species_pool, n_local - length(specialists))
  communities_conservation[i, c(specialists, generalists)] <- 1
}

# Fragmentos medianos (sitios 5-8): diversidad moderada
for(i in 5:8) {
  n_local <- sample(8:12, 1)
  generalists <- sample(11:n_species_pool, n_local)
  communities_conservation[i, generalists] <- 1
}

# Fragmentos pequeños (sitios 9-12): baja diversidad, solo generalistas
for(i in 9:12) {
  n_local <- sample(4:8, 1)
  generalists <- sample(15:n_species_pool, n_local)
  communities_conservation[i, generalists] <- 1
}

rownames(communities_conservation) <- paste("Fragment", 1:n_fragments)
colnames(communities_conservation) <- paste("Sp", 1:n_species_pool)

# Análisis de diversidad beta para conservación
beta_conservation <- beta.multi(communities_conservation)

print("Diversidad beta en paisaje fragmentado:")
## [1] "Diversidad beta en paisaje fragmentado:"
print(beta_conservation)
## $beta.SIM
## [1] 0.7002967
## 
## $beta.SNE
## [1] 0.09273605
## 
## $beta.SOR
## [1] 0.7930328
# Calcular contribución de cada sitio a la diversidad beta
contrib_sites <- numeric(n_fragments)
for(i in 1:n_fragments) {
  # Beta diversidad sin el sitio i
  comm_without_i <- communities_conservation[-i, ]
  beta_without_i <- beta.multi(comm_without_i)
  
  # Contribución = diferencia en beta total
  contrib_sites[i] <- beta_conservation$beta.SOR - beta_without_i$beta.SOR
}

# Crear ranking de prioridad
priority_ranking <- data.frame(
  Fragment = rownames(communities_conservation),
  Richness = rowSums(communities_conservation),
  Beta_contribution = contrib_sites,
  Priority_score = scale(rowSums(communities_conservation)) + scale(contrib_sites)
) %>%
  arrange(desc(Priority_score))

# Visualizar prioridades
ggplot(priority_ranking, aes(x = Richness, y = Beta_contribution, 
                            size = Priority_score, label = Fragment)) +
  geom_point(alpha = 0.7, color = "darkblue") +
  geom_text(hjust = -0.1, vjust = 0.5, size = 3) +
  theme_minimal() +
  labs(x = "Riqueza de especies", y = "Contribución a diversidad beta",
       title = "Priorización de fragmentos para conservación",
       size = "Puntaje de\nprioridad") +
  scale_size_continuous(range = c(2, 8))

print("Ranking de prioridad para conservación:")
## [1] "Ranking de prioridad para conservación:"
kable(priority_ranking, digits = 3)
Fragment Richness Beta_contribution Priority_score
Fragment 4 Fragment 4 14 0.016 1.905
Fragment 3 Fragment 3 14 0.013 1.190
Fragment 1 Fragment 1 16 0.008 0.420
Fragment 2 Fragment 2 17 0.007 0.309
Fragment 5 Fragment 5 8 0.016 0.216
Fragment 11 Fragment 11 5 0.018 0.172
Fragment 9 Fragment 9 7 0.015 -0.241
Fragment 12 Fragment 12 7 0.014 -0.373
Fragment 7 Fragment 7 9 0.012 -0.399
Fragment 10 Fragment 10 8 0.012 -0.713
Fragment 8 Fragment 8 11 0.007 -1.222
Fragment 6 Fragment 6 10 0.008 -1.264

11 Consideraciones Metodológicas

11.1 Limitaciones y precauciones

11.1.1 1. Efecto del esfuerzo de muestreo

# Demostrar efecto del esfuerzo de muestreo
set.seed(789)

# Comunidad verdadera
true_community <- rbinom(20, 1, 0.6)  # 20 especies, prob 0.6 cada una

# Simular diferentes esfuerzos de detección
detection_probs <- c(0.3, 0.5, 0.7, 0.9)
n_replicates <- 50

results_sampling <- data.frame()

for(det_prob in detection_probs) {
  beta_values <- numeric(n_replicates)
  
  for(rep in 1:n_replicates) {
    # Dos muestreos con diferente esfuerzo
    sample1 <- true_community * rbinom(20, 1, det_prob)
    sample2 <- true_community * rbinom(20, 1, det_prob)
    
    # Calcular beta diversidad
    if(sum(sample1) > 0 & sum(sample2) > 0) {
      beta_val <- vegdist(rbind(sample1, sample2), method = "jaccard")
      beta_values[rep] <- as.numeric(beta_val)
    } else {
      beta_values[rep] <- NA
    }
  }
  
  results_sampling <- rbind(results_sampling, 
                           data.frame(Detection_prob = det_prob,
                                     Beta_mean = mean(beta_values, na.rm = TRUE),
                                     Beta_sd = sd(beta_values, na.rm = TRUE)))
}

ggplot(results_sampling, aes(x = Detection_prob, y = Beta_mean)) +
  geom_line(size = 1.2, color = "red") +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = Beta_mean - Beta_sd, ymax = Beta_mean + Beta_sd),
                width = 0.02) +
  theme_minimal() +
  labs(x = "Probabilidad de detección", y = "Diversidad Beta promedio",
       title = "Efecto del esfuerzo de muestreo en diversidad beta",
       subtitle = "Comparación de muestras de la misma comunidad verdadera")

11.1.2 2. Selección de índices apropiados

# Comparación de índices bajo diferentes escenarios
scenarios <- list(
  "Recambio puro" = list(
    site1 = c(1,1,1,0,0,0),
    site2 = c(0,0,0,1,1,1)
  ),
  "Anidamiento puro" = list(
    site1 = c(1,1,1,1,1,1),
    site2 = c(1,1,1,0,0,0)
  ),
  "Mixto" = list(
    site1 = c(1,1,1,1,0,0),
    site2 = c(1,0,0,1,1,1)
  )
)

# Inicializar lista para almacenar resultados
results_list <- list()

for(scenario_name in names(scenarios)) {
  scenario <- scenarios[[scenario_name]]
  comm_matrix <- rbind(scenario$site1, scenario$site2)
  
  # SOLUCIÓN: Manejo seguro de cálculos de índices
  tryCatch({
    # Diferentes índices básicos
    jaccard <- vegdist(comm_matrix, method = "jaccard")
    sorensen <- vegdist(comm_matrix, method = "bray")  # equivalente para binarios
    
    # Inicializar valores por defecto
    jaccard_val <- ifelse(length(jaccard) > 0, as.numeric(jaccard), NA)
    sorensen_val <- ifelse(length(sorensen) > 0, as.numeric(sorensen), NA)
    turnover_val <- NA
    nestedness_val <- NA
    
    # Intentar descomposición de Baselga
    tryCatch({
      beta_decomp <- beta.pair(comm_matrix, index.family = "sorensen")
      
      # Verificar componentes disponibles y extraer valores de forma segura
      if("beta.sim" %in% names(beta_decomp) && !is.null(beta_decomp$beta.sim) && length(beta_decomp$beta.sim) > 0) {
        turnover_val <- as.numeric(beta_decomp$beta.sim)
      }
      
      if("beta.nes" %in% names(beta_decomp) && !is.null(beta_decomp$beta.nes) && length(beta_decomp$beta.nes) > 0) {
        nestedness_val <- as.numeric(beta_decomp$beta.nes)
      }
      
    }, error = function(e) {
      if(scenario_name == names(scenarios)[1]) {  # Solo mostrar mensaje una vez
        print(paste("Advertencia: No se pudo calcular descomposición de Baselga:", e$message))
      }
    })
    
    # Crear data.frame para este escenario
    scenario_result <- data.frame(
      Scenario = scenario_name,
      Jaccard = jaccard_val,
      Sorensen = sorensen_val,
      Turnover = turnover_val,
      Nestedness = nestedness_val,
      stringsAsFactors = FALSE
    )
    
    # Agregar a la lista de resultados
    results_list[[scenario_name]] <- scenario_result
    
  }, error = function(e) {
    print(paste("Error procesando escenario", scenario_name, ":", e$message))
    
    # Crear resultado con NAs en caso de error
    results_list[[scenario_name]] <<- data.frame(
      Scenario = scenario_name,
      Jaccard = NA,
      Sorensen = NA,
      Turnover = NA,
      Nestedness = NA,
      stringsAsFactors = FALSE
    )
  })
}

# SOLUCIÓN: Combinar resultados de forma segura
if(length(results_list) > 0) {
  # Verificar que todos los data.frames tienen la misma estructura
  all_same_structure <- all(sapply(results_list, function(x) {
    identical(names(x), names(results_list[[1]]))
  }))
  
  if(all_same_structure) {
    # Combinar usando do.call para mayor seguridad
    comparison_indices <- do.call(rbind, results_list)
    rownames(comparison_indices) <- NULL  # Limpiar nombres de filas
  } else {
    print("Error: Estructura inconsistente entre escenarios")
    # Crear tabla básica con estructura estándar
    comparison_indices <- data.frame(
      Scenario = names(scenarios),
      Jaccard = rep(NA, length(scenarios)),
      Sorensen = rep(NA, length(scenarios)),
      Turnover = rep(NA, length(scenarios)),
      Nestedness = rep(NA, length(scenarios))
    )
  }
} else {
  print("Error: No se pudieron procesar ninguno de los escenarios")
  comparison_indices <- data.frame(
    Scenario = character(0),
    Jaccard = numeric(0),
    Sorensen = numeric(0),
    Turnover = numeric(0),
    Nestedness = numeric(0)
  )
}

# Mostrar tabla de resultados
if(nrow(comparison_indices) > 0) {
  kable(comparison_indices, digits = 3, 
        caption = "Comparación de índices bajo diferentes escenarios")
  
  # Visualizar diferencias solo si hay datos válidos
  # Filtrar columnas con datos válidos (no todos NA)
  cols_to_plot <- names(comparison_indices)[-1]  # Excluir columna Scenario
  valid_cols <- cols_to_plot[sapply(cols_to_plot, function(col) {
    !all(is.na(comparison_indices[[col]]))
  })]
  
  if(length(valid_cols) > 0) {
    # Transformar a formato largo para ggplot
    comparison_long <- comparison_indices %>%
      select(Scenario, all_of(valid_cols)) %>%
      pivot_longer(cols = -Scenario, names_to = "Index", values_to = "Value") %>%
      filter(!is.na(Value))

    if(nrow(comparison_long) > 0) {
      ggplot(comparison_long, aes(x = Scenario, y = Value, fill = Index)) +
        geom_col(position = "dodge") +
        scale_fill_brewer(palette = "Set2") +
        theme_minimal() +
        labs(y = "Valor del índice", 
             title = "Comportamiento de diferentes índices",
             subtitle = "Comparación bajo escenarios contrastantes") +
        theme(axis.text.x = element_text(angle = 45, hjust = 1))
    } else {
      print("No hay datos válidos para visualizar")
      plot.new()
      text(0.5, 0.5, "No hay datos válidos para la visualización", cex = 1.2)
    }
  } else {
    print("No se pudieron calcular índices válidos")
    plot.new()
    text(0.5, 0.5, "Error: No se pudieron calcular índices válidos", cex = 1.2)
  }
} else {
  print("No hay resultados para mostrar")
}

# Diagnóstico adicional
print("Diagnóstico de cálculos:")
## [1] "Diagnóstico de cálculos:"
print(paste("Número de escenarios procesados:", length(results_list)))
## [1] "Número de escenarios procesados: 3"
print(paste("Filas en tabla final:", ifelse(exists("comparison_indices"), nrow(comparison_indices), 0)))
## [1] "Filas en tabla final: 3"

12 Conclusiones y Perspectivas Futuras

12.1 Síntesis de conceptos clave

La diversidad beta representa uno de los componentes más informativos del análisis de biodiversidad, proporcionando insights cruciales sobre:

  1. Procesos ecológicos: Ensamblaje de comunidades, dispersión, especiación
  2. Patrones espaciales: Efectos de la fragmentación, gradientes ambientales
  3. Conservación: Priorización de áreas, diseño de redes de protección
  4. Monitoreo: Detección de cambios temporales en biodiversidad

12.2 Avances metodológicos recientes

12.2.1 Análisis funcional y filogenético

# Ejemplo conceptual de diversidad beta funcional
# library(FD)
# library(ape)

# Calcular diversidad beta basada en rasgos funcionales
# functional_beta <- beta.pair(comm_matrix, trait_matrix)

# Incorporar información filogenética
# phylo_beta <- beta.pair.abund(comm_matrix, phylo_tree)

12.2.2 Modelos jerárquicos Bayesianos

La integración de enfoques Bayesianos permite: - Incorporar incertidumbre en estimaciones - Modelar efectos jerárquicos complejos - Propagación de errores de muestreo

12.2.3 Big data y sensores remotos

  • eDNA: Detección de especies mediante ADN ambiental
  • Sensores acústicos: Monitoreo continuo de comunidades
  • Imágenes satelitales: Estimación de diversidad a gran escala

12.3 Recomendaciones para investigación futura

  1. Estandarización de protocolos de muestreo y análisis
  2. Integración de múltiples dimensiones (taxonómica, funcional, filogenética)
  3. Desarrollo de métricas específicas para datos temporales
  4. Incorporación de efectos de cambio climático

12.4 Referencias principales

  • Baselga, A. (2010). Partitioning the turnover and nestedness components of beta diversity. Global Ecology and Biogeography, 19(1), 134-143.

  • Legendre, P., & De Cáceres, M. (2013). Beta diversity as the variance of community data: dissimilarity coefficients and partitioning. Ecology Letters, 16(8), 951-963.

  • Whittaker, R. H. (1960). Vegetation of the Siskiyou mountains, Oregon and California. Ecological Monographs, 30(3), 279-338.

  • Anderson, M. J., et al. (2011). Navigating the multiple meanings of β diversity: a roadmap for the practicing ecologist. Ecology Letters, 14(1), 19-28.

Documento generado con R Markdown. Para más información sobre los paquetes utilizados, consulte la documentación de CRAN.