Código
install.packages(c("tidyverse", "vegan", "iNEXT", "betapart", "car", "ggrepel"))Inteligencia Artificial aplicada a los datos en Ciencias Biológicas
| Módulo | 3 de 6 — Estadística biológica: fundamentos y multivariada |
| Intensidad del módulo | 3.5 horas (curso completo: 16 horas) |
| Modalidad | Virtual asincrónica (autoguiada), con sesiones en vivo de acompañamiento |
| Dirigido a | Profesionales egresados en Biología |
| Herramientas | R, RStudio, tidyverse, vegan, iNEXT, betapart y asistentes de IA |
| Institución | Universidad de La Guajira — Facultad de Ciencias Básicas y Aplicadas, Programa de Biología |
| Coordinación académica | Dr. Cristian Enrique Granados Martínez |
Sobre el autor
Cristian Enrique Granados Martínez es biólogo y Doctor en Ciencias Naturales para el Desarrollo (DOCINADE, Instituto Tecnológico de Costa Rica), con énfasis en Gestión de Recursos Naturales. Es docente-investigador del Programa de Biología de la Universidad de La Guajira, donde coordina el Comité de Investigación y Extensión de la Facultad de Ciencias Básicas y Aplicadas, e integra los grupos de investigación en Biotecnología y en Ecología y Biodiversidad de Ecosistemas Tropicales (EBET). Es Investigador Asociado reconocido por Minciencias.
Con más de 17 años de experiencia de campo en la Sierra Nevada de Santa Marta y La Guajira, su investigación se centra en los macroinvertebrados acuáticos, la limnología y la biodiversidad de agua dulce del Caribe colombiano. Ha descrito cuatro especies nuevas para la ciencia y es autor de más de 28 artículos indexados, 4 libros y 17 capítulos de libro.
ORCID: 0000-0003-0161-9537 · Sitio web: biolocristiam.github.io
Cómo citar esta guía: Granados-Martínez, C. E. (2026). Guía de estadística biológica con IA como copiloto. Módulo 3: Estadística biológica, fundamentos y multivariada. Curso Inteligencia Artificial aplicada a los datos en Ciencias Biológicas. Universidad de La Guajira, Riohacha, Colombia.
1 Antes de empezar
1.1 Cómo usar esta guía
En el Módulo 2 convertiste una libreta de campo en dos objetos listos para analizar: datos_integrados y matriz_comunidad. En este módulo les haces preguntas biológicas. ¿Difieren las estaciones? ¿Qué gradiente ordena las comunidades? ¿Qué dice la fauna sobre la calidad del agua? ¿Cuánta diversidad se pierde, y de qué tipo?
El énfasis está en entender cada método: qué pregunta responde, qué supone y cómo se lee su resultado. La IA te ayudará a escribir e interpretar, pero solo puedes corregirla si sabes qué debería responder. En el Módulo 4 convertirás lo aprendido aquí en recetas asistidas por IA.
| Bloque | Contenido | Tiempo |
|---|---|---|
| 3.1 | Estadística descriptiva e inferencial | ≈ 45 min |
| 3.2 | Multivariada: PCA, clúster, NMDS, RDA y CCA | ≈ 75 min |
| 3.3 | Fundamentos de bioindicación | ≈ 45 min |
| 3.4 | Fundamentos de análisis de diversidad | ≈ 45 min |
Encontrarás los mismos recuadros del Módulo 2:
Tip💡 Consejo
Buenas prácticas y atajos que te ahorran tiempo.
Importante🤖 Prompt para tu copiloto
Una instrucción lista para copiar en tu asistente de IA. Úsala, modifícala y compara lo que te responde con lo que dice esta guía.
Advertencia⚠️ Cuidado
Errores frecuentes o situaciones en las que la IA suele equivocarse.
Nota✅ Punto de control
Ejercicio para verificar lo aprendido, con solución desplegable.
1.2 Objetivos de aprendizaje
Al terminar este módulo podrás:
- Describir variables ecológicas y explicar por qué los conteos rara vez son normales.
- Elegir una prueba de hipótesis, verificar sus supuestos e interpretar su resultado.
- Ejecutar e interpretar PCA, clúster, NMDS, RDA y CCA.
- Explicar cómo se construyen BMWP/Col, ASPT y EPT, y cuáles son sus limitaciones.
- Comparar comunidades con números de Hill, rarefacción y extrapolación.
- Separar la diversidad beta en recambio y anidamiento.
1.3 Paquetes
Instala una sola vez los paquetes nuevos de este módulo:
Y cárgalos al inicio de tu documento:
Código
library(tidyverse) # manejo de datos y gráficos
library(vegan) # ordenación y análisis de comunidades
library(iNEXT) # rarefacción y extrapolación
library(betapart) # partición de la diversidad beta
library(car) # prueba de Levene
library(ggrepel) # etiquetas que no se superponen
theme_set(theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold"),
legend.position = "bottom"))
colores_epoca <- c("Seca" = "#D9822B", "Lluvias" = "#2B6CB0")1.4 Los datos del curso
Seguimos con el mismo caso: seis estaciones de un río de la vertiente norte de la Sierra Nevada de Santa Marta, muestreadas en época seca y de lluvias, con tres réplicas (36 muestras). E1 es la cabecera, a 870 m, y E6 la desembocadura, a 15 m. La estación E4 está aguas abajo de una descarga municipal.
Advertencia⚠️ Datos simulados
Los datos son simulados con fines didácticos. Reproducen patrones realistas, pero no corresponden a un muestreo real y no deben citarse como tales. Las interpretaciones de esta guía valen solo como ejercicio.
Para que la guía funcione por sí sola, el siguiente bloque reconstruye los objetos del Módulo 2, ya limpios. Puedes desplegarlo si quieres verlo, pero no necesitas estudiarlo.
Ver el código que reconstruye los datos del Módulo 2
set.seed(2026)
estaciones <- tibble(
estacion = paste0("E", 1:6),
nombre = c("Cabecera", "Bosque ripario", "Finca cafetera",
"Bajo descarga municipal", "Recuperación", "Desembocadura"),
altitud_m = c(870, 610, 395, 210, 95, 15),
latitud = c(11.0312, 11.0687, 11.1045, 11.1418, 11.1897, 11.2364),
longitud = c(-73.4521, -73.4703, -73.4852, -73.4987, -73.5125, -73.5281),
uso_suelo = c("Bosque", "Bosque", "Agrícola", "Urbano", "Agrícola", "Mixto")
)
perturbacion <- c(E1 = 0.05, E2 = 0.15, E3 = 0.35, E4 = 0.85, E5 = 0.55, E6 = 0.45)
muestras <- expand_grid(estacion = estaciones$estacion,
epoca = c("Seca", "Lluvias"), replica = 1:3) |>
left_join(estaciones |> select(estacion, altitud_m), by = "estacion") |>
mutate(perturbacion = perturbacion[estacion],
muestra = sprintf("%s-%s-R%d", estacion,
ifelse(epoca == "Seca", "S", "L"), replica))
taxa <- tribble(
~orden, ~familia, ~sens, ~base,
"Ephemeroptera", "Baetidae", 0.60, 60,
"Ephemeroptera", "Leptohyphidae", 0.65, 35,
"Ephemeroptera", "Leptophlebiidae", 0.85, 30,
"Ephemeroptera", "Caenidae", 0.40, 12,
"Plecoptera", "Perlidae", 0.95, 14,
"Trichoptera", "Hydropsychidae", 0.60, 40,
"Trichoptera", "Philopotamidae", 0.85, 10,
"Trichoptera", "Hydrobiosidae", 0.85, 6,
"Trichoptera", "Glossosomatidae", 0.80, 8,
"Trichoptera", "Leptoceridae", 0.70, 9,
"Trichoptera", "Calamoceratidae", 0.90, 5,
"Trichoptera", "Hydroptilidae", 0.55, 8,
"Coleoptera", "Elmidae", 0.70, 30,
"Coleoptera", "Psephenidae", 0.90, 6,
"Coleoptera", "Ptilodactylidae", 0.85, 5,
"Odonata", "Libellulidae", 0.40, 7,
"Odonata", "Gomphidae", 0.70, 5,
"Odonata", "Coenagrionidae", 0.35, 8,
"Megaloptera", "Corydalidae", 0.75, 4,
"Hemiptera", "Naucoridae", 0.60, 7,
"Hemiptera", "Veliidae", 0.55, 8,
"Diptera", "Simuliidae", 0.55, 35,
"Diptera", "Chironomidae", 0.15, 120,
"Diptera", "Tipulidae", 0.45, 5,
"Diptera", "Ceratopogonidae", 0.30, 8,
"Diptera", "Psychodidae", 0.10, 12,
"Diptera", "Culicidae", 0.05, 9,
"Basommatophora", "Physidae", 0.10, 25,
"Basommatophora", "Planorbidae", 0.25, 8,
"Amphipoda", "Hyalellidae", 0.50, 10,
"Decapoda", "Palaemonidae", 0.55, 6,
"Haplotaxida", "Tubificidae", 0.02, 50,
"Rhynchobdellida", "Glossiphoniidae", 0.15, 6
)
orden_estaciones <- paste0("E", 1:6)
macro <- expand_grid(muestras, taxa) |>
mutate(tolerancia = plogis(12 * ((1 - sens) - perturbacion)),
enriquec = ifelse(sens < 0.3, 0.25 + 2.5 * perturbacion, 1),
lambda = base * tolerancia * enriquec *
ifelse(epoca == "Lluvias", 0.55, 1),
abundancia = rnbinom(n(), mu = lambda, size = 1.2)) |>
filter(abundancia > 0) |>
select(muestra, estacion, epoca, replica, orden, familia, abundancia) |>
mutate(estacion = factor(estacion, levels = orden_estaciones),
epoca = factor(epoca, levels = c("Seca", "Lluvias")))
fq <- muestras |>
mutate(
temperatura_c = round(19.5 + (870 - altitud_m) / 870 * 8.5 +
ifelse(epoca == "Seca", 0.8, -0.6) + rnorm(n(), 0, 0.5), 1),
od_mg_l = round(8.4 - 5.2 * perturbacion +
ifelse(epoca == "Lluvias", 0.4, 0) + rnorm(n(), 0, 0.35), 2),
ph = round(7.1 + 0.6 * perturbacion + rnorm(n(), 0, 0.15), 2),
conductividad_us_cm = round(45 + 420 * perturbacion^1.5 +
(870 - altitud_m) * 0.08 + rnorm(n(), 0, 12), 0),
turbidez_ntu = round(exp(log(3) + 2.2 * perturbacion +
ifelse(epoca == "Lluvias", 1.1, 0) + rnorm(n(), 0, 0.25)), 1),
velocidad_m_s = round(pmax(0.05, 0.85 - (870 - altitud_m) / 870 * 0.45 +
ifelse(epoca == "Lluvias", 0.25, 0) + rnorm(n(), 0, 0.08)), 2)
)
fq$od_mg_l[14] <- NA # sensor sin batería
fq$turbidez_ntu[27] <- NA
metricas_muestra <- macro |>
group_by(muestra, estacion, epoca) |>
summarise(
riqueza = n_distinct(familia),
abundancia_total = sum(abundancia),
riqueza_ept = n_distinct(familia[orden %in% c("Ephemeroptera", "Plecoptera", "Trichoptera")]),
prop_chironomidae = sum(abundancia[familia == "Chironomidae"]) / sum(abundancia),
.groups = "drop"
)
datos_integrados <- metricas_muestra |>
left_join(fq |> select(muestra, temperatura_c, od_mg_l, ph,
conductividad_us_cm, turbidez_ntu, velocidad_m_s),
by = "muestra") |>
left_join(estaciones |> mutate(estacion = factor(estacion, levels = orden_estaciones)),
by = "estacion")
matriz_comunidad <- macro |>
select(muestra, familia, abundancia) |>
pivot_wider(names_from = familia, values_from = abundancia, values_fill = 0) |>
arrange(muestra) |>
column_to_rownames("muestra")
matriz_estacion <- macro |>
group_by(estacion, familia) |>
summarise(abundancia = sum(abundancia), .groups = "drop") |>
pivot_wider(names_from = familia, values_from = abundancia, values_fill = 0) |>
column_to_rownames("estacion")Comprueba que tienes lo necesario:
Código
dim(datos_integrados) # 36 muestras[1] 36 18Código
dim(matriz_comunidad) # 36 muestras × 33 familias[1] 36 33Código
# ¿Las filas de ambas tablas están en el mismo orden?
identical(rownames(matriz_comunidad), datos_integrados$muestra)[1] TRUE
Advertencia⚠️ La verificación que casi nadie hace
vegan relaciona la comunidad y el ambiente por posición de fila, no por nombre. Si las tablas están en distinto orden, R no da error: calcula un resultado falso que parece correcto. La última línea debe devolver TRUE. Repítela cada vez que filtres o reordenes una tabla. La IA casi nunca incluye esta comprobación en su código.
2 Bloque 3.1 — Estadística descriptiva e inferencial
2.1 Describir antes de probar
Las medidas de tendencia central (media, mediana) dicen dónde está el valor típico. Las de dispersión (desviación estándar, rango intercuartílico, coeficiente de variación) dicen cuánto varían los datos alrededor de él.
Código
datos_integrados |>
group_by(estacion) |>
summarise(
riqueza_media = mean(riqueza),
riqueza_de = sd(riqueza),
abund_media = mean(abundancia_total),
abund_mediana = median(abundancia_total),
abund_cv_pct = 100 * sd(abundancia_total) / mean(abundancia_total)
) |>
mutate(across(where(is.numeric), \(x) round(x, 1)))Observa E4: su riqueza es menos de la mitad que en el resto, pero su abundancia no se desploma. La descarga cambió quién vive allí más que cuántos individuos hay.
2.2 Por qué los datos biológicos rara vez son normales
Los conteos por taxón tienen tres características que chocan con la distribución normal: no pueden ser negativos, tienen muchos ceros y unos pocos valores muy altos.
Código
conteos <- as.vector(as.matrix(matriz_comunidad))
pct_ceros <- 100 * mean(conteos == 0)
tibble(conteo = conteos) |>
ggplot(aes(conteo)) +
geom_histogram(bins = 40, fill = "#1F4E5F") +
labs(x = "Individuos por familia y muestra", y = "Frecuencia")
El 35 % de las celdas de la matriz son ceros, y la media de los conteos (9.9) queda muy por debajo de la varianza (808). Cuando la varianza supera ampliamente la media hablamos de sobredispersión, algo típico de organismos que se distribuyen en parches.
Tip💡 Las métricas agregadas se comportan mejor
Un conteo por familia es muy asimétrico, pero la riqueza o el índice BMWP de una muestra suelen estar mucho más cerca de la normalidad. Por eso los supuestos se revisan para cada variable y cada modelo. Si necesitas modelar conteos por taxón, la vía adecuada son los modelos lineales generalizados (GLM) con distribución binomial negativa.
2.3 Cómo elegir la prueba
| Pregunta | Datos cumplen supuestos | No los cumplen |
|---|---|---|
| ¿Difieren dos grupos? | t de Student (Welch) | Mann-Whitney |
| ¿Difieren tres grupos o más? | ANOVA + Tukey | Kruskal-Wallis + comparaciones por pares |
| ¿Se relacionan dos variables? | Pearson | Spearman |
Toda prueba parte de una hipótesis nula (H₀), normalmente “no hay diferencia”. El valor p es la probabilidad de obtener un resultado tan extremo como el observado si H₀ fuera cierta. No es la probabilidad de que H₀ sea cierta, y tampoco mide el tamaño ni la importancia biológica del efecto.
2.4 Dos grupos: t de Student y Mann-Whitney
¿El oxígeno disuelto es menor en E4 que en la cabecera?
Código
e1_e4 <- datos_integrados |>
filter(estacion %in% c("E1", "E4")) |>
mutate(estacion = droplevels(estacion))
t.test(od_mg_l ~ estacion, data = e1_e4) # Welch, por defecto
Welch Two Sample t-test
data: od_mg_l by estacion
t = 15.669, df = 6.9055, p-value = 1.188e-06
alternative hypothesis: true difference in means between group E1 and group E4 is not equal to 0
95 percent confidence interval:
3.274440 4.442226
sample estimates:
mean in group E1 mean in group E4
8.218333 4.360000 R aplica por defecto la corrección de Welch, que no exige varianzas iguales. La diferencia es de unos 3.9 mg/L, y ese número importa más que el valor p.
La alternativa no paramétrica compara rangos en lugar de valores. La usamos para otra pregunta: ¿difiere la proporción de Chironomidae entre las dos estaciones? Al ser una proporción acotada entre 0 y 1, conviene no suponer normalidad:
Código
wilcox.test(prop_chironomidae ~ estacion, data = e1_e4)
Wilcoxon rank sum exact test
data: prop_chironomidae by estacion
W = 3, p-value = 0.01515
alternative hypothesis: true location shift is not equal to 0
Advertencia⚠️ Con pocas réplicas, la prueba no paramétrica no puede ser significativa
Con 3 réplicas por grupo, el valor p más pequeño posible en Mann-Whitney es 0.10, incluso con los grupos totalmente separados. Aquí comparamos seis muestras por estación (dos épocas × tres réplicas), y el mínimo posible es 0.002. Antes de muestrear, calcula si tu diseño puede detectar un efecto.
2.5 Tres grupos o más: ANOVA
¿La riqueza depende de la estación, de la época o de su interacción?
Código
modelo_riqueza <- aov(riqueza ~ estacion * epoca, data = datos_integrados)
tabla_anova <- summary(modelo_riqueza)
tabla_anova Df Sum Sq Mean Sq F value Pr(>F)
estacion 5 1544.1 308.83 99.266 3.15e-15 ***
epoca 1 2.3 2.25 0.723 0.403
estacion:epoca 5 11.9 2.38 0.766 0.583
Residuals 24 74.7 3.11
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La estación tiene un efecto fuerte. La época y la interacción, no. Calculamos también qué parte de la variación explica cada término (η²):
Código
ss <- tabla_anova[[1]][["Sum Sq"]]
round(setNames(ss / sum(ss), trimws(rownames(tabla_anova[[1]]))), 3) estacion epoca estacion:epoca Residuals
0.946 0.001 0.007 0.046 2.5.1 Verificar los supuestos
ANOVA supone que los residuos son aproximadamente normales, con varianzas similares entre grupos, y que las observaciones son independientes.
Código
residuos <- residuals(modelo_riqueza)
shapiro.test(residuos) # normalidad
Shapiro-Wilk normality test
data: residuos
W = 0.98006, p-value = 0.7476Código
leveneTest(riqueza ~ estacion * epoca, data = datos_integrados) # homocedasticidadGráfico cuantil-cuantil de los residuos del ANOVA. Los puntos deben seguir la línea.
Código
ggplot(tibble(residuos), aes(sample = residuos)) +
stat_qq(color = "#1F4E5F") + stat_qq_line() +
labs(x = "Cuantiles teóricos", y = "Residuos")
Ninguna de las dos pruebas rechaza su hipótesis nula y el gráfico no muestra desviaciones graves, así que el ANOVA es razonable.
Advertencia⚠️ Dos errores que la IA comete a menudo
- Probar la normalidad de los datos crudos. Lo que debe ser normal son los residuos. Si mezclas seis estaciones muy distintas, los datos crudos pueden verse “no normales” aunque el modelo sea perfecto.
- Tratar un valor p > 0.05 como prueba de normalidad. Con pocas réplicas, Shapiro-Wilk tiene poca potencia. Mira siempre el gráfico.
2.5.2 Comparaciones entre pares
El ANOVA dice que alguna estación difiere, no cuál. La prueba de Tukey compara todos los pares y controla el error por hacer muchas comparaciones. Estos son los pares que no difieren:
Código
TukeyHSD(modelo_riqueza, which = "estacion")$estacion |>
as.data.frame() |>
rownames_to_column("par") |>
filter(`p adj` > 0.05) |>
mutate(across(where(is.numeric), \(x) round(x, 3)))Las estaciones forman niveles de riqueza a lo largo del río, y E4 no se parece a ninguna otra.
2.6 Cuando no se cumplen los supuestos: Kruskal-Wallis
La abundancia total varía mucho más en unas estaciones que en otras:
Código
leveneTest(abundancia_total ~ estacion, data = datos_integrados)Levene rechaza la igualdad de varianzas, de modo que usamos Kruskal-Wallis:
Código
kruskal.test(abundancia_total ~ estacion, data = datos_integrados)
Kruskal-Wallis rank sum test
data: abundancia_total by estacion
Kruskal-Wallis chi-squared = 4.2896, df = 5, p-value = 0.5085No hay evidencia de diferencias en abundancia total, aunque la riqueza cambia drásticamente. Es una lección ecológica: la abundancia total, por sí sola, no detecta esta perturbación, porque los organismos tolerantes reemplazan a los sensibles.
Advertencia⚠️ ¿Qué podemos concluir? Pseudorreplicación
Las tres réplicas de una estación vienen del mismo tramo y del mismo día. Permiten comparar estas seis estaciones, pero no son réplicas independientes del factor “contaminación”. Tenemos un río y una descarga, así que no podemos generalizar a “las descargas municipales” (Hurlbert, 1984). Un copiloto escribirá esa conclusión general si no lo detienes.
Importante🤖 Prompt para tu copiloto
Hice un ANOVA de dos vías (estación × época) con riqueza de familias:
6 estaciones de un solo río, 2 épocas y 3 réplicas por combinación.
[pegar tabla_anova y las pruebas de supuestos]
Redacta un párrafo de resultados para un artículo. Después dime qué
conclusiones NO permite este diseño y por qué.Revisa que informe el tamaño del efecto, no solo el valor p, y que mencione los límites del diseño.
Nota✅ Punto de control 3.1
Ejercicio. Compara la riqueza EPT (riqueza_ept) entre estaciones. Verifica los supuestos, elige la prueba adecuada y di qué pares difieren de E4.
Solución:
Código
modelo_ept <- aov(riqueza_ept ~ estacion, data = datos_integrados)
shapiro.test(residuals(modelo_ept))
Shapiro-Wilk normality test
data: residuals(modelo_ept)
W = 0.95532, p-value = 0.1538Código
leveneTest(riqueza_ept ~ estacion, data = datos_integrados)Código
summary(modelo_ept) Df Sum Sq Mean Sq F value Pr(>F)
estacion 5 401.1 80.23 41.62 1.31e-12 ***
Residuals 30 57.8 1.93
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Código
# Pares que incluyen a E4
TukeyHSD(modelo_ept)$estacion |>
as.data.frame() |>
rownames_to_column("par") |>
filter(str_detect(par, "E4")) |>
mutate(across(where(is.numeric), \(x) round(x, 3)))Código
# Comprobación no paramétrica
kruskal.test(riqueza_ept ~ estacion, data = datos_integrados)
Kruskal-Wallis rank sum test
data: riqueza_ept by estacion
Kruskal-Wallis chi-squared = 30.114, df = 5, p-value = 1.4e-05Código
pairwise.wilcox.test(datos_integrados$riqueza_ept, datos_integrados$estacion,
p.adjust.method = "holm", exact = FALSE)
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: datos_integrados$riqueza_ept and datos_integrados$estacion
E1 E2 E3 E4 E5
E2 0.294 - - - -
E3 0.066 0.066 - - -
E4 0.066 0.066 0.066 - -
E5 0.066 0.066 0.294 0.066 -
E6 0.066 0.066 0.868 0.066 0.243
P value adjustment method: holm Ni Shapiro ni Levene rechazan sus hipótesis, así que ANOVA y Tukey son adecuados, y E4 difiere de todas las estaciones. Kruskal-Wallis llega a la misma conclusión global. Observa las comparaciones por pares: con seis muestras por grupo y 15 comparaciones, la corrección de Holm no deja bajar el valor p ajustado de cerca de 0.07, aunque los grupos estén totalmente separados. Elegir una prueba no paramétrica “por seguridad” tiene un costo real en potencia.
3 Bloque 3.2 — Estadística multivariada: ordenación y clasificación
3.1 La lógica general
Con 33 familias y 36 muestras no podemos graficar cada variable por separado. Los métodos multivariados reducen muchas variables a unos pocos ejes o grupos interpretables.
| Método | Datos de entrada | Pregunta | Tipo |
|---|---|---|---|
| PCA | Variables ambientales | ¿Qué gradientes ambientales hay? | Ordenación indirecta |
| Clúster | Disimilitudes | ¿Qué muestras forman grupos? | Clasificación |
| NMDS | Disimilitudes de comunidades | ¿Qué tan parecidas son las comunidades? | Ordenación indirecta |
| RDA | Comunidad + ambiente | ¿Cuánta composición explica el ambiente? (respuesta lineal) | Ordenación directa |
| CCA | Comunidad + ambiente | Igual, con respuesta unimodal | Ordenación directa |
La ordenación indirecta busca primero el patrón de los datos y luego lo relaciona con otras variables. La directa solo busca la parte del patrón que explican las variables ambientales que le das.
3.2 PCA: los gradientes ambientales
El PCA crea nuevos ejes, combinaciones lineales de las variables originales, que recogen la mayor varianza posible. Antes de aplicarlo:
- Se eliminan las filas con
NA, porque el PCA no los admite. - Se transforman las variables muy asimétricas (log de conductividad y turbidez).
- Se estandarizan todas (
scale. = TRUE). Si no, la conductividad, en cientos de µS/cm, dominaría sobre el pH.
Código
ambiente <- datos_integrados |>
select(muestra, estacion, epoca, temperatura_c, od_mg_l, ph,
conductividad_us_cm, turbidez_ntu, velocidad_m_s) |>
drop_na() |>
mutate(log_cond = log10(conductividad_us_cm),
log_turb = log10(turbidez_ntu))
nrow(ambiente) # perdimos 2 muestras por datos faltantes[1] 34Código
vars_amb <- ambiente |>
select(temperatura_c, od_mg_l, ph, log_cond, log_turb, velocidad_m_s)
pca <- prcomp(vars_amb, scale. = TRUE)
summary(pca)$importance[, 1:3] |> round(3) # varianza explicada PC1 PC2 PC3
Standard deviation 2.114 1.015 0.458
Proportion of Variance 0.745 0.172 0.035
Cumulative Proportion 0.745 0.917 0.952Código
round(pca$rotation[, 1:2], 2) # cargas PC1 PC2
temperatura_c 0.43 0.28
od_mg_l -0.44 0.00
ph 0.42 -0.25
log_cond 0.47 -0.04
log_turb 0.33 -0.67
velocidad_m_s -0.35 -0.65Código
sitios_pca <- ambiente |>
select(muestra, estacion, epoca) |>
bind_cols(as_tibble(pca$x[, 1:2]))
flechas_pca <- as_tibble(pca$rotation[, 1:2] * 3.2, rownames = "variable")
var_pca <- 100 * summary(pca)$importance[2, 1:2]
ggplot(sitios_pca, aes(PC1, PC2)) +
geom_hline(yintercept = 0, color = "grey85") +
geom_vline(xintercept = 0, color = "grey85") +
geom_point(aes(color = estacion, shape = epoca), size = 3) +
geom_segment(data = flechas_pca, aes(x = 0, y = 0, xend = PC1, yend = PC2),
arrow = arrow(length = unit(0.2, "cm")), color = "grey30") +
geom_text_repel(data = flechas_pca, aes(label = variable), color = "grey20") +
scale_color_viridis_d(end = 0.9) +
labs(x = sprintf("PC1 (%.1f %%)", var_pca[1]),
y = sprintf("PC2 (%.1f %%)", var_pca[2]),
color = "Estación", shape = "Época")
Cómo leerlo. Dos ejes retienen el 92 % de la varianza. PC1 enfrenta oxígeno con conductividad y temperatura: separa las aguas frías, oxigenadas y poco mineralizadas de la cabecera de las estaciones bajas y alteradas. PC2 está cargado por turbidez y velocidad, y separa la época seca de la de lluvias.
Advertencia⚠️ Un eje no es una causa
En este río, altitud y perturbación avanzan casi juntas. PC1 mezcla ambos gradientes, y el PCA no puede separarlos. Nombrar un eje “contaminación” es una interpretación tuya, que debes justificar con conocimiento del sistema.
3.3 Clúster jerárquico: agrupar comunidades
El clúster aplica dos decisiones, y ambas afectan al resultado:
- Medida de disimilitud. Bray-Curtis es el estándar para abundancias: va de 0 (comunidades idénticas) a 1 (nada en común) e ignora las ausencias compartidas. Eso es deseable, porque que dos sitios carezcan de una familia no los hace parecidos.
- Método de agrupamiento. Usamos el promedio (UPGMA). El método de Ward se suele pedir por costumbre, pero supone distancias euclidianas, y Bray-Curtis no lo es.
La raíz cuadrada reduce el peso de las familias dominantes, como Chironomidae:
Código
dist_bray <- vegdist(sqrt(matriz_comunidad), method = "bray")
cluster <- hclust(dist_bray, method = "average")
plot(cluster, hang = -1, cex = 0.7, main = "", xlab = "", sub = "",
ylab = "Disimilitud de Bray-Curtis")
rect.hclust(cluster, k = 3, border = c("#1F4E5F", "#D9822B", "#2B6CB0"))
Código
grupos <- cutree(cluster, k = 3)
table(grupo = grupos, estacion = datos_integrados$estacion) estacion
grupo E1 E2 E3 E4 E5 E6
1 6 6 1 0 0 0
2 0 0 5 0 6 6
3 0 0 0 6 0 0Código
coef_cofenetico <- cor(dist_bray, cophenetic(cluster))
coef_cofenetico[1] 0.7447537Aparecen tres comunidades: las estaciones de arriba (E1 y E2, junto con una muestra de E3), las intermedias (E3, E5, E6) y E4 aparte. El coeficiente cofenético (0.74) mide qué tan bien el dendrograma conserva las disimilitudes originales. Se suele considerar aceptable por encima de 0.7.
3.4 NMDS: ordenar comunidades
El NMDS coloca las muestras en un plano de modo que las más parecidas queden cerca. Trabaja con el orden de las disimilitudes, no con sus valores, y por eso exige muy poco de los datos. Los ejes no tienen varianza explicada y su rotación es arbitraria.
Código
set.seed(2026) # el NMDS parte de configuraciones aleatorias
nmds <- metaMDS(sqrt(matriz_comunidad), distance = "bray", k = 2,
trymax = 100, autotransform = FALSE, trace = FALSE)
nmds$stress[1] 0.1322014El stress mide cuánto se distorsionan las disimilitudes al llevarlas a dos dimensiones:
| Stress | Calidad de la representación (Clarke, 1993) |
|---|---|
| < 0.05 | Excelente |
| < 0.10 | Buena |
| < 0.20 | Utilizable, con cautela |
| > 0.30 | Casi arbitraria |
Nuestro stress de 0.132 indica una representación útil. Añadimos las variables ambientales con envfit(), que las proyecta sobre la ordenación ya hecha:
Código
filas_amb <- match(ambiente$muestra, rownames(matriz_comunidad))
set.seed(2026)
ajuste_amb <- envfit(scores(nmds, display = "sites")[filas_amb, ],
vars_amb, permutations = 999)
sitios_nmds <- as_tibble(scores(nmds, display = "sites"), rownames = "muestra") |>
left_join(datos_integrados |> select(muestra, estacion, epoca), by = "muestra")
poligonos <- sitios_nmds |>
group_by(estacion) |>
slice(chull(NMDS1, NMDS2)) |>
ungroup()
flechas_nmds <- as_tibble(scores(ajuste_amb, display = "vectors") * 0.65,
rownames = "variable") |>
filter(ajuste_amb$vectors$pvals < 0.05)
ggplot(sitios_nmds, aes(NMDS1, NMDS2)) +
geom_polygon(data = poligonos, aes(fill = estacion), alpha = 0.18) +
geom_point(aes(color = estacion, shape = epoca), size = 3) +
geom_segment(data = flechas_nmds, aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
arrow = arrow(length = unit(0.2, "cm")), color = "grey30") +
geom_text_repel(data = flechas_nmds, aes(label = variable), color = "grey20") +
scale_color_viridis_d(end = 0.9) +
scale_fill_viridis_d(end = 0.9) +
labs(color = "Estación", fill = "Estación", shape = "Época") +
coord_equal()
Las comunidades se ordenan en un gradiente, con E1 y E2 en un extremo y E4 en el otro. Las estaciones E5 y E6 no vuelven a la posición de la cabecera: la recuperación aguas abajo es parcial.
3.4.1 ¿Las diferencias son reales? PERMANOVA
El gráfico sugiere grupos, pero necesitamos una prueba. PERMANOVA compara la disimilitud dentro y entre grupos, y obtiene el valor p por permutaciones, sin suponer normalidad (Anderson, 2001).
Código
set.seed(2026)
permanova <- adonis2(dist_bray ~ estacion + epoca, data = datos_integrados,
permutations = 999, by = "margin")
permanovaLa estación explica el 55 % de la variación en composición, y la época el 4 %. Hay un supuesto clave: PERMANOVA es sensible a las diferencias de dispersión. Comprobamos que los grupos tengan una dispersión similar (Anderson, 2006):
Código
set.seed(2026)
permutest(betadisper(dist_bray, datos_integrados$estacion), permutations = 999)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 5 0.015960 0.0031920 1.0884 999 0.397
Residuals 30 0.087981 0.0029327 Un resultado no significativo indica que la diferencia detectada se debe a la posición de los grupos y no a su dispersión.
3.5 RDA y CCA: cuánto explica el ambiente
3.5.1 ¿Lineal o unimodal?
La RDA supone que las especies responden de forma lineal al gradiente y la CCA que tienen un óptimo (respuesta unimodal). Una regla práctica es mirar la longitud del primer eje de un análisis de correspondencia sin tendencia (DCA):
Código
dca <- decorana(matriz_comunidad)
largo_dca <- diff(range(scores(dca, display = "sites", choices = 1)))
round(largo_dca, 2)[1] 2.23Como guía, menos de 3 unidades favorece RDA y más de 4 favorece CCA. Nuestro gradiente mide unas 2.2 unidades, así que usamos RDA. Es una orientación, no una ley.
Los conteos crudos no son adecuados para RDA, porque usa distancias euclidianas. La transformación de Hellinger lo resuelve y evita tratar las ausencias compartidas como parecido (Legendre y Gallagher, 2001).
3.5.2 Colinealidad
Si incluyes variables redundantes, los coeficientes se vuelven inestables. El factor de inflación de la varianza (VIF) lo detecta, y por encima de 10 suele indicar un problema.
Código
Y <- decostand(matriz_comunidad[ambiente$muestra, ], method = "hellinger")
identical(rownames(Y), ambiente$muestra) # comprobación obligatoria[1] TRUECódigo
rda_completo <- rda(Y ~ od_mg_l + log_cond + log_turb + temperatura_c +
velocidad_m_s + ph, data = ambiente)
round(vif.cca(rda_completo), 1) od_mg_l log_cond log_turb temperatura_c velocidad_m_s
6.4 22.9 5.0 10.7 7.3
ph
4.5 La conductividad y la temperatura superan el umbral. Retiramos la conductividad, la variable más redundante, y también el pH, muy correlacionado con el oxígeno y la turbidez, para tener un modelo más sencillo:
Código
rda_final <- rda(Y ~ od_mg_l + log_turb + temperatura_c + velocidad_m_s,
data = ambiente)
round(vif.cca(rda_final), 1) od_mg_l log_turb temperatura_c velocidad_m_s
3.7 3.2 7.5 7.0
Advertencia⚠️ Quitar una variable no la descarta como causa
Retirar la conductividad no significa que no importe. Significa que su información se solapa con la del oxígeno y que este modelo no puede separarlas. La elección es tuya y debe justificarse biológicamente, no delegarse en un criterio automático.
3.5.3 Significancia e interpretación
Código
r2_rda <- RsquareAdj(rda_final)
round(unlist(r2_rda), 3) r.squared adj.r.squared
0.470 0.397 Código
set.seed(2026)
anova(rda_final, permutations = 999) # modelo globalCódigo
set.seed(2026)
anova(rda_final, by = "margin", permutations = 999) # efecto único de cada variableEl ambiente explica el 40 % de la composición (R² ajustado). El oxígeno es el único con efecto marginal significativo. Eso no quiere decir que las demás variables no importen: la prueba marginal mide lo que cada variable aporta después de las otras, y aquí comparten gran parte de la información.
Código
sitios_rda <- as_tibble(scores(rda_final, display = "wa", choices = 1:2, scaling = 2),
rownames = "muestra") |>
left_join(ambiente |> select(muestra, estacion, epoca), by = "muestra")
familias_rda <- as_tibble(scores(rda_final, display = "species", choices = 1:2, scaling = 2),
rownames = "familia") |>
mutate(peso = sqrt(RDA1^2 + RDA2^2)) |>
slice_max(peso, n = 10)
flechas_rda <- as_tibble(scores(rda_final, display = "bp", choices = 1:2, scaling = 2),
rownames = "variable")
escala <- 0.8 * max(abs(sitios_rda$RDA1)) / max(abs(flechas_rda$RDA1))
flechas_rda <- flechas_rda |> mutate(across(c(RDA1, RDA2), \(x) x * escala))
var_rda <- 100 * eigenvals(rda_final)[1:2] / rda_final$tot.chi
ggplot(sitios_rda, aes(RDA1, RDA2)) +
geom_hline(yintercept = 0, color = "grey85") +
geom_vline(xintercept = 0, color = "grey85") +
geom_point(aes(color = estacion, shape = epoca), size = 3) +
geom_segment(data = flechas_rda, aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
arrow = arrow(length = unit(0.2, "cm")), color = "grey30") +
geom_text_repel(data = flechas_rda, aes(label = variable), color = "grey20") +
geom_text_repel(data = familias_rda, aes(label = familia), size = 3.2,
fontface = "italic", color = "#8A4B14", max.overlaps = Inf) +
scale_color_viridis_d(end = 0.9) +
labs(x = sprintf("RDA1 (%.1f %% de la variación total)", var_rda[1]),
y = sprintf("RDA2 (%.1f %% de la variación total)", var_rda[2]),
color = "Estación", shape = "Época")
Cómo leer el triplot (escalamiento 2).
- Una familia apuntando en la dirección de una flecha tiende a aumentar con esa variable.
- El ángulo entre dos flechas refleja su correlación: agudo si es positiva, recto si es nula, obtuso si es negativa.
- Las muestras proyectadas sobre una flecha indican su valor aproximado en esa variable.
Aquí, Tubificidae y Chironomidae aumentan hacia el extremo de E4, opuesto al oxígeno, mientras familias como Baetidae, Leptohyphidae e Hydropsychidae se asocian a las estaciones oxigenadas. El primer eje concentra casi toda la variación explicada.
3.5.4 CCA
Si tu gradiente fuera largo, el código cambia poco. La CCA usa la matriz de conteos, sin Hellinger:
Código
cca_final <- cca(matriz_comunidad[ambiente$muestra, ] ~
od_mg_l + log_turb + temperatura_c + velocidad_m_s,
data = ambiente)
set.seed(2026)
anova(cca_final, permutations = 999)
Advertencia⚠️ Cuidado con las familias raras en la CCA
La CCA se basa en la distancia chi-cuadrado, que da mucho peso a los taxones raros. Una sola familia en una sola muestra puede arrastrar un eje.
Importante🤖 Prompt para tu copiloto
Te adjunto la tabla de VIF, el R² ajustado, las pruebas anova() global
y marginal, y las coordenadas de una RDA con transformación de Hellinger
de macroinvertebrados acuáticos [pegar resultados].
Interpreta qué gradiente representa RDA1. Separa con claridad lo que
muestran los datos de lo que es hipótesis tuya, y dime qué variable
de confusión no puedo descartar con este diseño.Comprueba que no confunda correlación con causalidad y que mencione la relación entre altitud y perturbación.
Nota✅ Punto de control 3.2
Ejercicio. Repite el NMDS y el PERMANOVA usando solo presencia/ausencia, con la disimilitud de Jaccard. ¿Cambia el patrón? ¿Por qué podría cambiar?
Solución:
Código
dist_jaccard <- vegdist(matriz_comunidad, method = "jaccard", binary = TRUE)
set.seed(2026)
nmds_pa <- metaMDS(dist_jaccard, k = 2, trymax = 100, trace = FALSE)
nmds_pa$stress[1] 0.09306454Código
set.seed(2026)
adonis2(dist_jaccard ~ estacion + epoca, data = datos_integrados,
permutations = 999, by = "margin")NMDS con presencia/ausencia (Jaccard).
Código
as_tibble(scores(nmds_pa, display = "sites"), rownames = "muestra") |>
left_join(datos_integrados |> select(muestra, estacion, epoca), by = "muestra") |>
ggplot(aes(NMDS1, NMDS2, color = estacion, shape = epoca)) +
geom_point(size = 3) +
scale_color_viridis_d(end = 0.9) +
labs(color = "Estación", shape = "Época") +
coord_equal()
La separación entre estaciones se mantiene porque E4 pierde muchas familias. La época, en cambio, deja de ser significativa: su efecto está sobre todo en cuántos individuos hay, no en qué familias aparecen. Con presencia/ausencia, la pregunta cambia de “¿cuántos?” a “¿quiénes?”.
4 Bloque 3.3 — Fundamentos de bioindicación
4.1 Qué es un bioindicador
Un bioindicador es un organismo o conjunto de organismos cuya presencia, abundancia o condición refleja el estado del ambiente. La fisicoquímica es una fotografía del agua en el momento del muestreo. La fauna integra las condiciones de semanas o meses, y conserva la señal de un vertido puntual que quizá ya pasó.
Los macroinvertebrados acuáticos son un caso clásico por varias razones:
- Viven en el agua durante todo o gran parte de su ciclo de vida.
- Se mueven poco, así que reflejan las condiciones locales.
- Incluyen grupos con tolerancias muy distintas.
- Son fáciles de muestrear y, a nivel de familia, de identificar.
4.2 Cómo se construye un índice biótico
4.2.1 BMWP/Col
El BMWP nació en Reino Unido (Armitage et al., 1983) y Roldán (2003) lo adaptó a Colombia. La lógica es sencilla:
- Cada familia recibe un puntaje de 1 a 10 según su tolerancia a la contaminación orgánica. Un 10 es muy sensible y un 1, muy tolerante.
- Se identifican las familias presentes en la muestra.
- Se suman sus puntajes. La abundancia no interviene.
Código
puntajes_bmwp <- tribble(
~familia, ~puntaje,
"Calamoceratidae", 10,
"Gomphidae", 10,
"Perlidae", 10,
"Psephenidae", 10,
"Ptilodactylidae", 10,
"Hydrobiosidae", 9,
"Leptophlebiidae", 9,
"Philopotamidae", 9,
"Leptoceridae", 8,
"Palaemonidae", 8,
"Simuliidae", 8,
"Veliidae", 8,
"Baetidae", 7,
"Caenidae", 7,
"Coenagrionidae", 7,
"Glossosomatidae", 7,
"Hyalellidae", 7,
"Hydropsychidae", 7,
"Hydroptilidae", 7,
"Leptohyphidae", 7,
"Naucoridae", 7,
"Psychodidae", 7,
"Corydalidae", 6,
"Elmidae", 6,
"Libellulidae", 6,
"Planorbidae", 5,
"Ceratopogonidae", 3,
"Glossiphoniidae", 3,
"Physidae", 3,
"Tipulidae", 3,
"Chironomidae", 2,
"Culicidae", 2,
"Tubificidae", 1
)
# ¿Todas las familias de nuestros datos tienen puntaje?
setdiff(unique(macro$familia), puntajes_bmwp$familia)character(0)
Advertencia⚠️ Verifica la tabla en la fuente original
Estos puntajes siguen la tabla de Roldán (2003), reproducida en Roldán-Pérez (2016). La tabla publicada no es única. Álvarez (2005) presenta una versión modificada, por ejemplo con Psychodidae en 2 y Simuliidae en 5. La reproducción de 2016 incluye a Hydropsychidae con 7 y con 5; aquí usamos 7. Una IA te dará puntajes con total seguridad, a veces mezclando versiones. Cita siempre qué tabla usaste.
Calculamos el índice por muestra. El ASPT divide el BMWP entre el número de familias con puntaje y estima la sensibilidad promedio, por lo que depende menos de la riqueza.
Código
bmwp_muestra <- macro |>
left_join(puntajes_bmwp, by = "familia") |>
group_by(muestra, estacion, epoca) |>
summarise(bmwp = sum(puntaje),
familias = n_distinct(familia),
aspt = bmwp / familias,
.groups = "drop")
head(bmwp_muestra)El valor se traduce a una clase de calidad. Esta es la tabla que se suele citar de Roldán (2003). Observa que deja un hueco entre 121 y 150, así que aquí tratamos todo valor superior a 100 como clase “Buena”. Las versiones publicadas también difieren en estos límites, por lo que conviene indicar cuál se usó:
| Clase | Calidad | BMWP/Col | Significado |
|---|---|---|---|
| I | Buena | > 150; 101–120 | Aguas muy limpias o no contaminadas |
| II | Aceptable | 61–100 | Efectos evidentes de contaminación |
| III | Dudosa | 36–60 | Aguas moderadamente contaminadas |
| IV | Crítica | 16–35 | Aguas muy contaminadas |
| V | Muy crítica | < 15 | Aguas fuertemente contaminadas |
Código
clasificar_bmwp <- function(x) {
cut(x, breaks = c(-Inf, 15, 35, 60, 100, Inf),
labels = c("Muy crítica", "Crítica", "Dudosa", "Aceptable", "Buena"))
}
bmwp_muestra <- bmwp_muestra |>
mutate(clase = clasificar_bmwp(bmwp))
count(bmwp_muestra, estacion, clase)4.2.2 EPT
Los órdenes Ephemeroptera, Plecoptera y Trichoptera agrupan muchas familias sensibles. Se suele reportar su riqueza o su proporción de individuos:
Código
ept_muestra <- macro |>
group_by(muestra) |>
summarise(pct_ept = 100 * sum(abundancia[orden %in% c("Ephemeroptera", "Plecoptera", "Trichoptera")]) /
sum(abundancia))
indices <- bmwp_muestra |>
left_join(ept_muestra, by = "muestra") |>
left_join(datos_integrados |>
select(muestra, riqueza_ept, od_mg_l, conductividad_us_cm),
by = "muestra")
indices |>
group_by(estacion) |>
summarise(across(c(bmwp, aspt, riqueza_ept, pct_ept), mean)) |>
mutate(across(where(is.numeric), \(x) round(x, 1)))Código
ggplot(indices, aes(estacion, bmwp)) +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 100, ymax = Inf, fill = "#2B6CB0", alpha = 0.10) +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 60, ymax = 100, fill = "#4C9F70", alpha = 0.10) +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 35, ymax = 60, fill = "#E6B800", alpha = 0.12) +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = 15, ymax = 35, fill = "#D9822B", alpha = 0.14) +
annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 15, fill = "#B03A2E", alpha = 0.14) +
geom_point(aes(color = epoca), size = 3, alpha = 0.8,
position = position_dodge(width = 0.5)) +
scale_color_manual(values = colores_epoca) +
labs(x = "Estación (cabecera → desembocadura)", y = "BMWP/Col", color = "Época")
4.3 Relación con la fisicoquímica
Un índice es útil si responde al gradiente de calidad del agua:
Código
indices |>
pivot_longer(c(od_mg_l, conductividad_us_cm), names_to = "variable", values_to = "valor") |>
mutate(variable = if_else(variable == "od_mg_l", "Oxígeno disuelto (mg/L)",
"Conductividad (µS/cm)")) |>
ggplot(aes(valor, bmwp)) +
geom_point(aes(color = estacion), size = 2.5) +
scale_color_viridis_d(end = 0.9) +
facet_wrap(~ variable, scales = "free_x") +
labs(x = NULL, y = "BMWP/Col", color = "Estación")
Código
cor.test(~ bmwp + od_mg_l, data = indices, method = "spearman", exact = FALSE)
Spearman's rank correlation rho
data: bmwp and od_mg_l
S = 594.62, p-value = 1.062e-14
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.9167192 Código
cor.test(~ bmwp + conductividad_us_cm, data = indices, method = "spearman", exact = FALSE)
Spearman's rank correlation rho
data: bmwp and conductividad_us_cm
S = 15007, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.9314495 La asociación es muy fuerte, pero observa la nube de puntos: las muestras se agrupan por estación, de modo que no son observaciones independientes del gradiente. La correlación describe el patrón, pero su valor p es optimista.
4.4 Del número a la conclusión
El índice no se interpreta solo. Mira qué pasa al calcularlo con el mismo dato, pero uniendo todas las muestras de cada estación:
Código
bmwp_estacion <- macro |>
distinct(estacion, familia) |>
left_join(puntajes_bmwp, by = "familia") |>
group_by(estacion) |>
summarise(bmwp_unido = sum(puntaje), familias = n(),
aspt_unido = round(bmwp_unido / familias, 2)) |>
left_join(indices |> group_by(estacion) |>
summarise(bmwp_medio = round(mean(bmwp), 1),
aspt_medio = round(mean(aspt), 2)),
by = "estacion") |>
mutate(clase_media = clasificar_bmwp(bmwp_medio),
clase_unido = clasificar_bmwp(bmwp_unido))
bmwp_estacionE4 tiene clase Crítica como promedio de sus muestras y Aceptable al unirlas. El agua es la misma: lo que cambió fue el esfuerzo, porque cada muestra añadida puede sumar una familia más. Algunos individuos aislados (una Baetidae, una Hydropsychidae) suman el mismo puntaje que una población establecida. En cambio, el ASPT se mueve mucho menos.
Advertencia⚠️ Limitaciones del BMWP
- Ignora la abundancia. Un individuo cuenta igual que mil. E4 suma puntos por visitantes ocasionales, quizá arrastrados.
- Depende del esfuerzo. Compara solo muestras con el mismo método y esfuerzo.
- Mide sobre todo contaminación orgánica. Puede no detectar metales, pesticidas o sedimentación.
- Resolución de familia. Dentro de una familia hay especies con tolerancias distintas. Chironomidae es el ejemplo clásico.
- Calibración regional. Una tabla ajustada en un lugar no se transfiere automáticamente a ríos de otra región o altitud. Una estación de planicie puede tener valores bajos sin estar contaminada.
Una conclusión defendible combina varias líneas de evidencia:
E4 presenta BMWP/Col y ASPT bajos, ausencia casi total de EPT y dominancia de Chironomidae y Tubificidae, junto con baja concentración de oxígeno y alta conductividad. El patrón es coherente con enriquecimiento orgánico asociado a la descarga municipal, aunque se necesitan muestreos adicionales para atribuirlo a esa causa.
Importante🤖 Prompt para tu copiloto
Para estos resultados por estación [pegar bmwp_estacion e indices resumidos],
redacta una interpretación de la calidad del agua. No inventes puntajes
ni categorías. Si necesitas una tabla de referencia, pregúntame cuál uso.
Señala qué resultados dependen del esfuerzo de muestreo.Si te da una tabla de puntajes, contrástala con la fuente. Es uno de los casos en que la IA más inventa, porque hay muchas versiones parecidas.
Nota✅ Punto de control 3.3
Ejercicio. ¿Qué tan robusta es la clasificación a la tabla usada? Cambia dos puntajes como en Álvarez (2005): Psychodidae de 7 a 2 y Simuliidae de 8 a 5. Recalcula BMWP y di cuántas muestras cambian de clase.
Solución:
Código
puntajes_alt <- puntajes_bmwp |>
mutate(puntaje = case_when(familia == "Psychodidae" ~ 2,
familia == "Simuliidae" ~ 5,
TRUE ~ puntaje))
bmwp_alt <- macro |>
left_join(puntajes_alt, by = "familia") |>
group_by(muestra, estacion) |>
summarise(bmwp_alt = sum(puntaje), .groups = "drop") |>
mutate(clase_alt = clasificar_bmwp(bmwp_alt))
bmwp_comparado <- indices |>
select(muestra, bmwp, clase) |>
left_join(bmwp_alt, by = "muestra")
bmwp_comparado |>
group_by(estacion) |>
summarise(bmwp_original = round(mean(bmwp), 1),
bmwp_cambiado = round(mean(bmwp_alt), 1),
cambian_de_clase = sum(clase != clase_alt))Con solo dos puntajes distintos, 2 de 36 muestras cambian de clase, todas cerca de un límite. Por eso conviene reportar el valor numérico además de la categoría y declarar la tabla empleada.
5 Bloque 3.4 — Fundamentos de análisis de diversidad
5.1 Por qué un solo número no basta
Compara estas dos comunidades:
| Comunidad | Familia A | Familia B | Familia C | Familia D |
|---|---|---|---|---|
| 1 | 25 | 25 | 25 | 25 |
| 2 | 97 | 1 | 1 | 1 |
Las dos tienen 4 familias y 100 individuos, pero no son igual de diversas. La riqueza no ve la dominancia. Además, los índices clásicos no se comparan de forma intuitiva: un Shannon de 3 no es “el doble de diverso” que uno de 1.5.
5.2 Números de Hill: un lenguaje común
Los números de Hill (Hill, 1973; Jost, 2006) expresan la diversidad como un número efectivo de especies: cuántas especies igualmente abundantes darían la misma diversidad. El orden q controla cuánto pesan las especies abundantes:
| Orden | Equivale a | Qué representa |
|---|---|---|
| q = 0 | Riqueza | Cuenta todas las especies por igual |
| q = 1 | Exponencial de Shannon | Especies “comunes” |
| q = 2 | Inverso de Simpson | Especies “dominantes” |
Ahora sí podemos decir “esta comunidad tiene el doble de diversidad efectiva”. Una comunidad con 4 especies perfectamente equitativas vale 4 en los tres órdenes, y cuanto más caen los valores al aumentar q, mayor es la dominancia.
Código
hill <- renyi(matriz_estacion, scales = c(0, 1, 2), hill = TRUE)
colnames(hill) <- c("q0", "q1", "q2")
round(hill, 1)Código
perfil <- renyi(matriz_estacion, scales = seq(0, 3, 0.25), hill = TRUE) |>
as.data.frame() |>
rownames_to_column("estacion") |>
pivot_longer(-estacion, names_to = "q", values_to = "diversidad") |>
mutate(q = as.numeric(q))
ggplot(perfil, aes(q, diversidad, color = estacion)) +
geom_line(linewidth = 1) +
scale_color_viridis_d(end = 0.9) +
labs(x = "Orden q", y = "Número efectivo de familias", color = "Estación")
E4 pierde diversidad en todos los órdenes y su perfil cae rápido: pasa de 15 familias observadas a unas 3 dominantes. Si dos curvas se cruzan, ninguna comunidad es más diversa en todos los sentidos y conviene decirlo.
5.3 Rarefacción y extrapolación con iNEXT
La riqueza crece con el número de individuos. Comparar valores crudos premia al sitio mejor muestreado, y en nuestros datos los totales por estación van de 1379 a 2526 individuos.
- La rarefacción interpola hacia abajo: estima la diversidad para una muestra menor.
- La extrapolación proyecta hacia arriba, con fiabilidad hasta aproximadamente el doble de la muestra observada (Chao et al., 2014).
- La cobertura de muestreo estima qué proporción de los individuos de la comunidad pertenece a especies ya registradas. Una cobertura de 0.99 indica que falta por detectar cerca del 1 %.
iNEXT (Hsieh et al., 2016) necesita las familias en filas y las estaciones en columnas, de modo que trasponemos la matriz:
Código
abundancias <- as.data.frame(t(matriz_estacion))
set.seed(2026)
curvas <- iNEXT(abundancias, q = c(0, 1, 2), datatype = "abundance",
endpoint = 2 * min(colSums(abundancias)), nboot = 50)
Tip💡 nboot
nboot = 50 acelera el cálculo durante el aprendizaje. Para publicar, usa 200 o más.
Código
curvas$iNextEst$size_based |>
mutate(q = paste0("q = ", Order.q),
tramo = if_else(Method == "Extrapolation", "Extrapolación", "Rarefacción")) |>
ggplot(aes(m, qD, color = Assemblage, fill = Assemblage)) +
geom_ribbon(aes(ymin = qD.LCL, ymax = qD.UCL), alpha = 0.15, color = NA) +
geom_line(aes(linetype = tramo), linewidth = 0.8) +
scale_color_viridis_d(end = 0.9) +
scale_fill_viridis_d(end = 0.9) +
scale_linetype_manual(values = c("Rarefacción" = "solid", "Extrapolación" = "dashed")) +
facet_wrap(~ q, scales = "free_y") +
labs(x = "Número de individuos", y = "Número efectivo de familias",
color = "Estación", fill = "Estación", linetype = NULL) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))
Tip💡 El atajo
ggiNEXT(curvas) produce un gráfico similar en una línea. Aquí usamos ggplot2 para controlar el diseño.
La comparación más justa se hace a igual cobertura, porque comparamos comunidades igualmente completas, no igual número de individuos (Chao y Jost, 2012):
Código
set.seed(2026)
cobertura <- estimateD(abundancias, q = c(0, 1, 2), datatype = "abundance",
base = "coverage", level = 0.99, nboot = 50)
cobertura |>
mutate(q = paste0("q = ", Order.q)) |>
ggplot(aes(Assemblage, qD, color = Assemblage)) +
geom_pointrange(aes(ymin = qD.LCL, ymax = qD.UCL)) +
scale_color_viridis_d(end = 0.9, guide = "none") +
facet_wrap(~ q, scales = "free_y") +
labs(x = "Estación", y = "Número efectivo de familias")
Si los intervalos de confianza del 95 % no se solapan, la diferencia es clara. El orden en todos los valores de q es E1 > E2 > E3 > E6 ≈ E5 > E4.
Advertencia⚠️ Familias, no especies
iNEXT y los números de Hill se pensaron para especies. A nivel de familia, las curvas se saturan rápido y la cobertura suele ser alta. Puedes comparar estaciones, pero no hables de “riqueza de especies” y ten presente que la resolución taxonómica oculta diversidad.
5.4 Diversidad beta: qué tan distintas son las comunidades
Whittaker distinguió tres componentes: alfa (la diversidad local), gamma (la regional) y beta (el cambio entre sitios). Una medida simple es cuántas “comunidades distintas” caben en la región:
Código
gamma_rio <- ncol(matriz_estacion)
alfa_media <- mean(specnumber(matriz_estacion))
beta_whittaker <- gamma_rio / alfa_media
round(c(gamma = gamma_rio, alfa_media = alfa_media, beta = beta_whittaker), 2) gamma alfa_media beta
33.00 27.83 1.19 El valor mínimo es 1, cuando todas las estaciones tienen exactamente las mismas familias. Nuestro 1.19 indica que, al unir todas las muestras de cada estación, las listas de familias se parecen bastante. Esa medida no distingue cómo difieren.
5.4.1 Recambio o anidamiento
Dos comunidades pueden diferir de dos maneras muy distintas (Baselga, 2010):
- Recambio (turnover): unas especies reemplazan a otras.
- Anidamiento (nestedness): una comunidad es un subconjunto empobrecido de la otra.
La distinción tiene implicaciones de manejo. Si hay recambio, cada sitio aporta especies propias y hay que proteger varios. Si hay anidamiento, los sitios pobres pierden especies y el sitio más rico concentra la diversidad. betapart separa estos componentes (Baselga y Orme, 2012):
Código
incidencia <- decostand(matriz_estacion, method = "pa")
beta_sor <- beta.pair(incidencia, index.family = "sorensen")
beta_ab <- beta.pair.abund(matriz_estacion, index.family = "bray")
beta.multi(incidencia, index.family = "sorensen") |> unlist() |> round(3)beta.SIM beta.SNE beta.SOR
0.000 0.304 0.304 Comparamos cada estación con la cabecera:
Código
extraer <- function(balance, gradiente, datos) {
tibble(estacion = rownames(as.matrix(balance))[-1],
`Recambio / balance` = as.matrix(balance)[-1, "E1"],
`Anidamiento / gradiente` = as.matrix(gradiente)[-1, "E1"],
datos = datos)
}
bind_rows(
extraer(beta_sor$beta.sim, beta_sor$beta.sne, "Presencia/ausencia (Sørensen)"),
extraer(beta_ab$beta.bray.bal, beta_ab$beta.bray.gra, "Abundancia (Bray-Curtis)")
) |>
pivot_longer(c(`Recambio / balance`, `Anidamiento / gradiente`),
names_to = "componente", values_to = "disimilitud") |>
ggplot(aes(estacion, disimilitud, fill = componente)) +
geom_col() +
scale_fill_manual(values = c("Recambio / balance" = "#1F4E5F",
"Anidamiento / gradiente" = "#D9822B")) +
facet_wrap(~ datos) +
labs(x = "Estación comparada con E1", y = "Disimilitud", fill = NULL)
Los dos paneles cuentan historias complementarias:
- Por presencia/ausencia, todo es anidamiento. E1 tiene las 33 familias del río, de modo que cada estación es, por construcción, un subconjunto suyo.
- Por abundancia, domina el balance. Las mismas familias cambian de peso: las sensibles se vuelven escasas y las tolerantes se multiplican.
Advertencia⚠️ Una familia con un individuo basta
Al unir 18 muestras por estación, un individuo aislado cuenta como presencia. Por eso la incidencia puede subestimar el recambio, y la abundancia suele dar una imagen más completa. Antes de interpretar, pregúntate qué significa una “presencia” en tus datos.
Importante🤖 Prompt para tu copiloto
Tengo estos resultados de diversidad por estación: números de Hill,
diversidad a cobertura 0.99 y partición de la beta en balance y gradiente
[pegar tablas]. Los datos están a nivel de FAMILIA.
Interpreta los resultados para un biólogo. Explica qué significa que
E4 tenga poca diversidad en q = 2 y qué implica para el manejo.
Avísame si alguna afirmación necesita datos a nivel de especie.Revisa que no hable de “especies” sin aclararlo, y que distinga el resultado por incidencia del resultado por abundancia.
Nota✅ Punto de control 3.4
Ejercicio. Compara la diversidad de la época seca y la época de lluvias. Calcula primero la riqueza observada y después la diversidad a igual cobertura. ¿El número de individuos influye en la comparación?
Solución:
Código
abund_epoca <- macro |>
group_by(epoca, familia) |>
summarise(abundancia = sum(abundancia), .groups = "drop") |>
pivot_wider(names_from = epoca, values_from = abundancia, values_fill = 0) |>
column_to_rownames("familia")
colSums(abund_epoca) # individuos por época Seca Lluvias
7603 4100 Código
colSums(abund_epoca > 0) # riqueza observada Seca Lluvias
33 33 Código
set.seed(2026)
estimateD(abund_epoca, q = c(0, 1, 2), datatype = "abundance",
base = "coverage", level = 0.99, nboot = 50) |>
mutate(across(c(m, SC, qD, qD.LCL, qD.UCL), \(x) round(x, 2)))La época seca tiene casi el doble de individuos, porque las lluvias arrastran organismos, pero la riqueza observada es la misma. Al estandarizar por cobertura, las diferencias entre épocas son pequeñas: el arrastre reduce la abundancia mucho más que la diversidad.
6 Cierre del módulo
6.1 Lo que aprendiste
6.2 ¿Qué método uso?
| Mi pregunta | Método |
|---|---|
| ¿Difiere una métrica entre grupos? | t, ANOVA o sus alternativas no paramétricas |
| ¿Qué gradientes ambientales hay? | PCA |
| ¿Qué muestras se agrupan? | Clúster |
| ¿Cómo se ordenan las comunidades? | NMDS |
| ¿Difiere la composición entre grupos? | PERMANOVA + betadisper() |
| ¿Cuánta composición explica el ambiente? | RDA (gradiente corto) o CCA (largo) |
| ¿Cuál es la calidad biológica del agua? | BMWP/Col, ASPT, EPT y otras evidencias |
| ¿Cuál comunidad es más diversa? | Números de Hill a igual cobertura |
| ¿Por qué difieren las comunidades? | Partición de la diversidad beta |
6.3 Lo que viene
En el Módulo 4 convertirás estos fundamentos en recetas con IA como copiloto: limpieza de matrices, cálculo de índices de bioindicación, análisis de diversidad con interpretación asistida y reportes reproducibles en Quarto. Como ya sabes qué debe salir, podrás detectar cuando el copiloto se equivoca.
6.4 Banco de prompts del módulo
| Situación | Prompt sugerido |
|---|---|
| Elegir una prueba | “Tengo esta variable [descripción], estos grupos [n por grupo] y este diseño [descripción]. ¿Qué prueba usarías, qué supuestos debo revisar y qué harías si no se cumplen?” |
| Revisar supuestos | “Estos son los resultados de Shapiro y Levene y el gráfico de residuos [pegar]. ¿Puedo usar ANOVA? No te limites al valor p.” |
| Elegir la ordenación | “Tengo una matriz de [n] muestras × [n] taxones y estas variables ambientales. Mi pregunta es [pregunta]. ¿PCA, NMDS, RDA o CCA? Justifica.” |
| Leer un biplot | “Estas son las coordenadas de sitios, especies y variables de mi RDA, escalamiento 2 [pegar]. Explica el gradiente y separa hechos de hipótesis.” |
| Revisar código | “Revisa este código de vegan y señala riesgos: orden de las filas, transformaciones, semillas y número de permutaciones: [código]” |
| Interpretar diversidad | “Interpreta estos números de Hill a igual cobertura [pegar]. Mis datos son a nivel de familia. Indica qué no se puede concluir.” |
6.5 Errores frecuentes y cómo leerlos
| Mensaje de error | Qué significa | Solución habitual |
|---|---|---|
infinite or missing values in 'x' |
prcomp() encontró NA |
Filtra con drop_na() y alinea la matriz con las muestras que quedan |
variable lengths differ |
La comunidad y el ambiente tienen distinto número de filas | Usa las mismas muestras en ambas tablas |
unused arguments al usar recode() |
car tiene su propia función recode(), que oculta la de dplyr |
Escribe dplyr::recode() |
Stress is (nearly) zero |
Datos insuficientes para el NMDS | Revisa si hay filas duplicadas o muy pocas muestras |
Best solution was not repeated |
metaMDS() no halló dos soluciones similares |
Aumenta trymax, por ejemplo a 200 |
you have empty rows |
Hay muestras sin individuos | Elimínalas de la matriz antes de analizar |
cannot compute exact p-value with ties |
Hay empates en una prueba de rangos | Añade exact = FALSE e infórmalo |
6.6 Referencias
- Álvarez, L. F. (2005). Metodología para la utilización de los macroinvertebrados acuáticos como indicadores de la calidad del agua. Instituto de Investigación de Recursos Biológicos Alexander von Humboldt.
- Anderson, M. J. (2001). A new method for non-parametric multivariate analysis of variance. Austral Ecology, 26, 32–46.
- Anderson, M. J. (2006). Distance-based tests for homogeneity of multivariate dispersions. Biometrics, 62, 245–253.
- Armitage, P. D., Moss, D., Wright, J. F. y Furse, M. T. (1983). The performance of a new biological water quality score system based on macroinvertebrates over a wide range of unpolluted running-water sites. Water Research, 17, 333–347.
- Baselga, A. (2010). Partitioning the turnover and nestedness components of beta diversity. Global Ecology and Biogeography, 19, 134–143.
- Baselga, A. y Orme, C. D. L. (2012). betapart: an R package for the study of beta diversity. Methods in Ecology and Evolution, 3, 808–812.
- Chao, A., Gotelli, N. J., Hsieh, T. C., Sander, E. L., Ma, K. H., Colwell, R. K. y Ellison, A. M. (2014). Rarefaction and extrapolation with Hill numbers. Ecological Monographs, 84, 45–67.
- Chao, A. y Jost, L. (2012). Coverage-based rarefaction and extrapolation. Ecology, 93, 2533–2547.
- Clarke, K. R. (1993). Non-parametric multivariate analyses of changes in community structure. Australian Journal of Ecology, 18, 117–143.
- Hill, M. O. (1973). Diversity and evenness: a unifying notation and its consequences. Ecology, 54, 427–432.
- Hsieh, T. C., Ma, K. H. y Chao, A. (2016). iNEXT: an R package for rarefaction and extrapolation of species diversity. Methods in Ecology and Evolution, 7, 1451–1456.
- Hurlbert, S. H. (1984). Pseudoreplication and the design of ecological field experiments. Ecological Monographs, 54, 187–211.
- Jost, L. (2006). Entropy and diversity. Oikos, 113, 363–375.
- Legendre, P. y Gallagher, E. D. (2001). Ecologically meaningful transformations for ordination of species data. Oecologia, 129, 271–280.
- Roldán, G. (2003). Bioindicación de la calidad del agua en Colombia: uso del método BMWP/Col. Editorial Universidad de Antioquia.
- Roldán-Pérez, G. (2016). Los macroinvertebrados como bioindicadores de la calidad del agua: cuatro décadas de desarrollo en Colombia y Latinoamérica. Revista de la Academia Colombiana de Ciencias Exactas, Físicas y Naturales, 40(155), 254–274. https://doi.org/10.18257/raccefyn.335
6.7 Información de la sesión
Código
sessionInfo()R version 4.3.3 (2024-02-29 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 26200)
Matrix products: default
locale:
[1] LC_COLLATE=Spanish_Spain.utf8 LC_CTYPE=Spanish_Spain.utf8
[3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C
[5] LC_TIME=Spanish_Spain.utf8
time zone: America/Bogota
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggrepel_0.9.6 car_3.1-2 carData_3.0-5 betapart_1.6
[5] iNEXT_3.0.1 vegan_2.6-10 lattice_0.22-6 permute_0.9-7
[9] lubridate_1.9.3 forcats_1.0.0 stringr_1.6.0 dplyr_1.2.0
[13] purrr_1.0.2 readr_2.1.5 tidyr_1.3.2 tibble_3.2.1
[17] ggplot2_4.0.2 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] fastmatch_1.1-6 gtable_0.3.6 xfun_0.52 htmlwidgets_1.6.4
[5] tzdb_0.4.0 vctrs_0.7.1 tools_4.3.3 doSNOW_1.0.20
[9] generics_0.1.4 parallel_4.3.3 cluster_2.1.8.2 pkgconfig_2.0.3
[13] Matrix_1.6-5 RColorBrewer_1.1-3 S7_0.2.1 lifecycle_1.0.5
[17] compiler_4.3.3 farver_2.1.2 minpack.lm_1.2-4 codetools_0.2-19
[21] htmltools_0.5.8.1 snow_0.4-4 yaml_2.3.10 pillar_1.11.0
[25] MASS_7.3-60.0.1 iterators_1.0.14 abind_1.4-5 foreach_1.5.2
[29] picante_1.8.2 nlme_3.1-164 tidyselect_1.2.1 digest_0.6.35
[33] stringi_1.8.3 reshape2_1.4.4 labeling_0.4.3 splines_4.3.3
[37] magic_1.6-1 fastmap_1.1.1 grid_4.3.3 cli_3.6.2
[41] magrittr_2.0.3 ape_5.8 withr_3.0.2 scales_1.4.0
[45] timechange_0.3.0 rmarkdown_2.29 hms_1.1.3 evaluate_1.0.5
[49] knitr_1.50 rcdd_1.6 viridisLite_0.4.2 geometry_0.5.0
[53] mgcv_1.9-1 itertools_0.1-3 rlang_1.1.7 Rcpp_1.0.12
[57] glue_1.7.0 rstudioapi_0.16.0 jsonlite_1.8.8 R6_2.6.1
[61] plyr_1.8.9