Este informe reúne los tres puntos del taller: el análisis de salinidad y biomasa del suelo, el estudio del consumo de oxígeno en moluscos, y la evaluación de biodiversidad de anfibios según el tipo de hábitat. A continuación se cargan todas las bibliotecas y bases de datos necesarias para los tres análisis.

require(table1)
require(ggplot2)
require(agricolae)
require(car)
require(broom)
require(knitr)
require(kableExtra)
library(dplyr)
require(plotly)
require(tidyr)
require(gt)
require(htmltools)

load("C:/Users/Juan Esteban/Downloads/Salinidad.RData")
load("~/Archivos R/moluscos.RData")
load("C:/Users/Juan Esteban/Downloads/Biodiversidad.RData")

Punto 1 - Datos Salinidad

a) Estadísticas descriptivas

cv <- sapply(Salinidad, sd) / sapply(Salinidad, mean) * 100

df_desc <- data.frame(
  Variable = names(Salinidad),
  Media    = round(sapply(Salinidad, mean), 2),
  DE       = round(sapply(Salinidad, sd), 2),
  CV       = round(cv, 1),
  Mediana  = round(sapply(Salinidad, median), 2),
  Minimo   = round(sapply(Salinidad, min), 2),
  Maximo   = round(sapply(Salinidad, max), 2)
)
colnames(df_desc) <- c("Variable", "Media", "DE", "CV (%)", "Mediana", "Mínimo", "Máximo")

knitr::kable(df_desc, row.names = FALSE, align = "c",
             caption = "Tabla 1. Resumen de estadísticas descriptivas por variable")
Tabla 1. Resumen de estadísticas descriptivas por variable
Variable Media DE CV (%) Mediana Mínimo Máximo
Biomasa 1082.17 546.29 50.5 991.83 369.82 2337.33
pH 4.61 1.25 27.2 4.45 3.20 7.45
Salinidad 30.27 3.72 12.3 30.00 24.00 38.00
Zinc 17.83 8.27 46.4 19.24 0.21 31.29
Potasio 797.38 297.58 37.3 773.30 350.73 1441.67

Biomasa La biomasa presentó una media de 1080 g y una desviación estándar de 546 g, con valores entre 370 y 2340 g. El CV fue de 50.5% indicando una variabilidad relativamente alta entre las muestras. Mediana fue de 992 g.

pH El pH presentó una media de 4.61, desviación estándar de 1.25 y valores entre 3.2 y 7.45. Su coeficiente de variación fue de 27.2%. Mediana de 4.45.

Salinidad La salinidad presentó una media de 30.3, desviación estándar de 3.72 y un rango entre 24 y 38. El CV fue de 12.3%, siendo la variable con menor variabilidad relativa. Mediana de 30.

Zinc Presentó una media de 17.8, desviación estándar de 8.27 y valores entre 0.2011 y 31.3. Su CV fue de 46.4%, indicando una alta variabilidad relativa. Mediana de 19.2.

Potasio El potasio presentó una media de 797, desviación estándar de 298 y valores entre 351 y 1440. El CV fue de 37.3%, por lo que presentó una variabilidad relativa intermedia-alta. Mediana de 773.

par(mfrow = c(2,3))
hist(Salinidad$Biomasa, main="Biomasa", xlab="gr")
hist(Salinidad$pH,        main="pH",        xlab="pH")
hist(Salinidad$Salinidad, main="Salinidad", xlab="Salinidad")
hist(Salinidad$Zinc,      main="Zinc",      xlab="Zinc")
hist(Salinidad$Potasio,   main="Potasio",   xlab="Potasio")

Biomasa Media mayor a la mediana. Distribución asimétrica a la derecha (positiva). Una moda principal en 500-1500 y un grupito aparte en 2000-2500.

pH Media ligeramente mayor a la mediana, aunque la presencia de valores elevados entre aproximadamente 7.1 y 7.45 incrementa la dispersión y genera cierta asimetría. Asimetría positiva. Varias modas: una en 3-3.5, otra en 4.5-5 y un grupo en 7+.

Salinidad Media prácticamente igual a la mediana. Casi simétrica. Varias modas: una en 24-26, otra en 28-30 y en 34-36 (bimodal).

Zinc Media menor a la mediana. Asimetría negativa. Una moda en 15-25, con un grupito aparte cerca de 0.

Potasio Media un poco mayor a la mediana. Asimetría positiva (leve). Varias modas: 400-600, 800-1000 y 1200-1400.

par(mfrow = c(2,3))
boxplot(Salinidad$Biomasa,main="Biomasa")
boxplot(Salinidad$pH,main="pH")
boxplot(Salinidad$Salinidad, main="Salinidad")
boxplot(Salinidad$Zinc,main="Zinc")
boxplot(Salinidad$Potasio, main="Potasio")
par(mfrow = c(1,1))

Biomasa La mediana está casi en el centro de la caja, pero el bigote de arriba es mucho más largo que el de abajo. Indica asimetría positiva. No hay valores atípicos. Las muestras más altas están en el extremo del bigote.

pH La caja va de 3.45 a 5.35 y la mediana (4.45) está ubicada en medio. El bigote inferior es corto y el superior largo. Asimetría positiva. No hay valores atípicos.

Salinidad La mediana (30) queda justo en el centro de la caja y los bigotes son parecidos, el largo siendo un poco más largo. Es casi simétrica. Sin atípicos.

Zinc Mediana más cerca de Q3 que de Q1, bigote de abajo termina en 9.4. Sumado al punto suelto, indica asimetría negativa. Valor atípico que representa a cinco muestras (0.21 a 0.37). Forman un grupo con características propias y no errores de medición.

Potasio Bigote superior mucho más largo que el inferior. Asimetría positiva leve. Dentro de la caja, mediana un poco más cerca de Q3, así que las señales no son del todo claras. Sin atípicos.

b) Relación entre variables

R2 = cor(Salinidad[, c("Biomasa", "pH", "Salinidad", "Zinc", "Potasio")])^2

knitr::kable(round(R2, 3), align = "c",
             caption = "Tabla 2. Matriz de R² entre las variables numéricas")
Tabla 2. Matriz de R² entre las variables numéricas
Biomasa pH Salinidad Zinc Potasio
Biomasa 1.000 0.861 0.004 0.611 0.005
pH 0.861 1.000 0.002 0.519 0.001
Salinidad 0.004 0.002 1.000 0.182 0.000
Zinc 0.611 0.519 0.182 1.000 0.006
Potasio 0.005 0.001 0.000 0.006 1.000

pH: explica casi toda la variación de la biomasa. Queda un 13.9% sin explicar, que es la dispersión de los puntos alrededor de la recta.

Zinc: también explica mucho, pero unos 25 puntos porcentuales menos que el pH.

Salinidad y potasio: explican menos del 1%. Con 45 muestras, un R² tan bajo no es significativo. No sirve para predecir la biomasa.

ggplot(Salinidad, aes(x=pH, y=Biomasa)) + geom_point() + geom_smooth(method="lm")

pH vs Biomasa Los puntos siguen la recta de cerca y el intervalo de confianza del 95% es angosto, es decir, recta bien estimada de poca incertidumbre.

A mayor pH, mayor biomasa.

El pH explica cerca del 86% de la biomasa (R²≈0.86), este dato se confirma más abajo.

ggplot(Salinidad, aes(x=Zinc, y=Biomasa)) + geom_point() + geom_smooth(method="lm")

Zinc vs Biomasa La relación es negativa: a más Zinc, menos biomasa. Su R²=0.61 es fuerte, pero más dispersa que la del pH.

Entre Zinc 10 y 31 los puntos bajan con más dispersión, y uno con Zinc ~20 y biomasa ~1890 destaca por encima.

Parte de esta relación puede ser indirecta: pH y Zinc están correlacionados entre sí (R²=-0.72), y son las mismas muestras en el extremo. No se puede decir que el Zinc reduzca la biomasa. Solo que ambas varían juntas.

ggplot(Salinidad, aes(x=Salinidad, y=Biomasa)) + geom_point() + geom_smooth(method="lm")

Salinidad vs Biomasa La recta es casi horizontal y el intervalo de confianza es muy ancho. No se distingue el cero.

Los puntos forman columnas verticales porque la salinidad solo toma valores enteros. En cada columna hay biomasas muy distintas: con salinidad 30, van de unos 370 a más de 2300 g.

La salinidad no explica la biomasa en este conjunto de datos.

ggplot(Salinidad, aes(x=Potasio, y=Biomasa)) + geom_point() + geom_smooth(method="lm")

Potasio vs Biomasa Tampoco hay tendencia: recta casi plana y banda ancha que incluye la horizontal.

Se ven dos zonas: 5 puntos arriba a la izquierda (Potasio ~450-570, biomasa mayor a 2100) y muchos abajo con biomasa baja y Potasio parecido. Con Potasio ~500, la biomasa puede ser de 370 a 2337 g, así que el Potasio no sirve para predecirla.

c) Categorización del pH y análisis de varianza

if ("pH_cat" %in% names(Salinidad)) Salinidad$pH_cat = NULL   
cortes = quantile(Salinidad$pH, probs = c(0, 1/3, 2/3, 1))
cortes
##        0% 33.33333% 66.66667%      100% 
##  3.200000  3.883333  4.900000  7.450000
pH_cat = cut(Salinidad$pH, breaks = cortes, include.lowest = TRUE,
             labels = c("Bajo", "Medio", "Alto"))
table(pH_cat)
## pH_cat
##  Bajo Medio  Alto 
##    15    15    15
df_tabla3 <- data.frame(
  `Categoría de pH` = c("Bajo", "Medio", "Alto"),
  Media = round(tapply(Salinidad$Biomasa, pH_cat, mean), 2),
  DE = round(tapply(Salinidad$Biomasa, pH_cat, sd), 2),
  Mediana = round(tapply(Salinidad$Biomasa, pH_cat, median), 2),
  Mínimo = round(tapply(Salinidad$Biomasa, pH_cat, min), 2),
  Máximo = round(tapply(Salinidad$Biomasa, pH_cat, max), 2)
)

knitr::kable(
  df_tabla3,
  row.names = FALSE,
  align = "c",
  caption = "Tabla 3. Biomasa según categoría de pH"
)
Tabla 3. Biomasa según categoría de pH
Categoría.de.pH Media DE Mediana Mínimo Máximo
Bajo 593.02 164.83 545.54 369.82 977.52
Medio 1048.11 244.91 1039.64 568.46 1491.28
Alto 1605.39 547.60 1422.84 765.28 2337.33

Biomasa según nivel de pH La biomasa media aumenta escalonadamente con el nivel de pH.

En Bajo y Medio la media y la mediana casi coinciden (distribuciones simétricas); en Alto la media (1610) supera bastante a la mediana (1420), señal de asimetría positiva, coherente con las 5 muestras extremas (biomasa > 2100) que tiene este grupo.

ggplot(Salinidad, aes(x=pH_cat, y=Biomasa, fill=pH_cat)) + geom_boxplot()

Boxplot por categoría de pH En el boxplot, las tres cajas casi no se traslapan y suben de izquierda a derecha. Bajo tiene un atípico leve (~978 g). La caja de Alto es mucho más ancha que las otras dos: va de ~1200 a ~2190, contra ~480-660 en Bajo y ~890-1200 en Medio. Esa diferencia de tamaño es la primera pista visual de que las varianzas no son iguales.

anova_sal = aov(Salinidad$Biomasa ~ pH_cat)

H0: las medias de biomasa son iguales en los tres niveles de pH.
H1: al menos un nivel difiere.

Como p<0.05, se rechaza H0: el nivel de pH genera diferencias significativas en la biomasa.

tab_anova <- summary(anova_sal)[[1]]

df_anova <- data.frame(
  Fuente  = trimws(rownames(tab_anova)),
  gl      = tab_anova$Df,
  SC      = round(tab_anova$`Sum Sq`, 2),
  CM      = round(tab_anova$`Mean Sq`, 2),
  F       = round(tab_anova$`F value`, 3),
  p_valor = signif(tab_anova$`Pr(>F)`, 4)
)
colnames(df_anova) <- c("Fuente", "gl", "SC", "CM", "F", "p-valor")

knitr::kable(df_anova, row.names = FALSE, align = "c",
             caption = "Tabla 4. ANOVA de una vía — Biomasa según categoría de pH")
Tabla 4. ANOVA de una vía — Biomasa según categoría de pH
Fuente gl SC CM F p-valor
pH_cat 2 7712683 3856341.6 29.893 0
Residuals 42 5418235 129005.6 NA NA 
df_dec_anova <- data.frame(
  Prueba      = "ANOVA de una vía",
  H0          = "Las medias de biomasa son iguales en los tres niveles de pH",
  H1          = "Al menos un nivel difiere",
  Estadistico = round(tab_anova$`F value`[1], 3),
  p_valor     = signif(tab_anova$`Pr(>F)`[1], 4),
  Decision    = ifelse(tab_anova$`Pr(>F)`[1] < 0.05, "Se rechaza H0", "No se rechaza H0")
)
colnames(df_dec_anova) <- c("Prueba", "H0", "H1", "Estadístico F", "p-valor", "Decisión")

knitr::kable(df_dec_anova, row.names = FALSE, align = "c",
             caption = "Tabla 5. Decisión — ANOVA de una vía")
Tabla 5. Decisión — ANOVA de una vía
Prueba H0 H1 Estadístico F p-valor Decisión
ANOVA de una vía Las medias de biomasa son iguales en los tres niveles de pH Al menos un nivel difiere 29.893 0 Se rechaza H0 
shapiro_sal <- shapiro.test(residuals(anova_sal))

df_shapiro <- data.frame(
  Prueba      = "Shapiro-Wilk",
  H0          = "Los residuales siguen una distribución normal",
  Estadistico = round(shapiro_sal$statistic, 4),
  p_valor     = signif(shapiro_sal$p.value, 4),
  Decision    = ifelse(shapiro_sal$p.value < 0.05, "Se rechaza H0", "No se rechaza H0")
)
colnames(df_shapiro) <- c("Prueba", "H0", "Estadístico W", "p-valor", "Decisión")

knitr::kable(df_shapiro, row.names = FALSE, align = "c",
             caption = "Tabla 6. Prueba de normalidad de los residuales (Shapiro-Wilk)")
Tabla 6. Prueba de normalidad de los residuales (Shapiro-Wilk)
Prueba H0 Estadístico W p-valor Decisión
Shapiro-Wilk Los residuales siguen una distribución normal 0.9732 0.3749 No se rechaza H0

H0: los residuales siguen una distribución normal.
Como p>0.05, no se rechaza H0: el supuesto se cumple.

qqnorm(residuals(anova_sal)); qqline(residuals(anova_sal))

Q-Q Plot El Q-Qplot lo confirma: la mayoría de los puntos, sobre todo en el centro, sigue de cerca la recta. Hay una desviación leve en ambas colas (algunos puntos por debajo de la recta a la izquierda y por encima a la derecha), lo que indica colas un poco más pesadas que una normal perfecta, pero no lo suficiente como para que Shapiro lo detecte como anormal. Es un resultado coherente y no contradictorio.

bartlett_sal <- bartlett.test(Salinidad$Biomasa ~ pH_cat)
levene_sal   <- leveneTest(Salinidad$Biomasa ~ pH_cat)

df_homog <- data.frame(
  Prueba      = c("Bartlett", "Levene"),
  H0          = rep("Las varianzas de biomasa son iguales en los tres niveles de pH", 2),
  Estadistico = c(round(bartlett_sal$statistic, 3), round(levene_sal$`F value`[1], 3)),
  gl          = c(paste(round(bartlett_sal$parameter, 2), collapse = ", "),
                   paste(levene_sal$Df, collapse = ", ")),
  p_valor     = c(signif(bartlett_sal$p.value, 4), signif(levene_sal$`Pr(>F)`[1], 4)),
  Decision    = c(ifelse(bartlett_sal$p.value < 0.05, "Se rechaza H0", "No se rechaza H0"),
                   ifelse(levene_sal$`Pr(>F)`[1] < 0.05, "Se rechaza H0", "No se rechaza H0"))
)
colnames(df_homog) <- c("Prueba", "H0", "Estadístico", "gl", "p-valor", "Decisión")

knitr::kable(df_homog, row.names = FALSE, align = "c",
             caption = "Tabla 7. Pruebas de homogeneidad de varianzas (Bartlett y Levene)")
Tabla 7. Pruebas de homogeneidad de varianzas (Bartlett y Levene)
Prueba H0 Estadístico gl p-valor Decisión
Bartlett Las varianzas de biomasa son iguales en los tres niveles de pH 20.084 2 4.35e-05 Se rechaza H0
Levene Las varianzas de biomasa son iguales en los tres niveles de pH 8.675 2, 42 7.02e-04 Se rechaza H0

H0: las tres varianzas son iguales.
Ambas pruebas dan p<0.05, así que se rechaza H0 en las dos: el supuesto no se cumple.

Que las dos pruebas coincidan le da más fuerza a la conclusión: no es un artefacto de una sola prueba.

La causa se ve en las DE: 165 (Bajo), 245 (Medio) y 548 (Alto). El grupo Alto es 3.3 veces más variable que Bajo, porque mezcla pH “normales” (5-5.6) con el grupo aparte de pH mayor a 7 (biomasa 2160-2337).

welch_sal <- oneway.test(Salinidad$Biomasa ~ pH_cat)

df_welch <- data.frame(
  Prueba      = "Welch (ANOVA robusto a varianzas desiguales)",
  H0          = "Las medias de biomasa son iguales en los tres niveles de pH",
  Estadistico = round(welch_sal$statistic, 3),
  gl_num      = round(welch_sal$parameter[1], 2),
  gl_den      = round(welch_sal$parameter[2], 2),
  p_valor     = signif(welch_sal$p.value, 4),
  Decision    = ifelse(welch_sal$p.value < 0.05, "Se rechaza H0", "No se rechaza H0")
)
colnames(df_welch) <- c("Prueba", "H0", "Estadístico F", "gl num.", "gl denom.", "p-valor", "Decisión")

knitr::kable(df_welch, row.names = FALSE, align = "c",
             caption = "Tabla 8. Prueba de Welch (ANOVA robusto)")
Tabla 8. Prueba de Welch (ANOVA robusto)
Prueba H0 Estadístico F gl num. gl denom. p-valor Decisión
Welch (ANOVA robusto a varianzas desiguales) Las medias de biomasa son iguales en los tres niveles de pH 34.727 2 24.74 1e-07 Se rechaza H0

El resultado confirma al ANOVA clásico que el pH afecta la biomasa.

lsd_sal <- LSD.test(anova_sal, "pH_cat", console = FALSE)

df_lsd <- data.frame(
  Grupo = rownames(lsd_sal$groups),
  Media = round(lsd_sal$groups[, 1], 2),
  Letra = lsd_sal$groups$groups
)

knitr::kable(df_lsd, row.names = FALSE, align = "c",
             caption = "Tabla 9. Prueba LSD — comparaciones múltiples de medias por nivel de pH")
Tabla 9. Prueba LSD — comparaciones múltiples de medias por nivel de pH
Grupo Media Letra
Alto 1605.39 a
Medio 1048.11 b
Bajo 593.02 c

Letras distintas para los tres grupos significa que los tres niveles de pH difieren significativamente entre sí.

Tabla resumen — Punto c)

resumen_c <- data.frame(
  Prueba = c("ANOVA (clásica)",
             "Shapiro-Wilk (normalidad)",
             "Bartlett (homocedasticidad)",
             "Levene (homocedasticidad)",
             "Welch (ANOVA robusto)"),
  Estadistico = c(
    round(tab_anova$`F value`[1], 3),
    round(shapiro_sal$statistic, 4),
    round(bartlett_sal$statistic, 3),
    round(levene_sal$`F value`[1], 3),
    round(welch_sal$statistic, 3)
  ),
  p_valor = c(
    signif(tab_anova$`Pr(>F)`[1], 4),
    signif(shapiro_sal$p.value, 4),
    signif(bartlett_sal$p.value, 4),
    signif(levene_sal$`Pr(>F)`[1], 4),
    signif(welch_sal$p.value, 4)
  ),
  Decision = c(
    ifelse(tab_anova$`Pr(>F)`[1] < 0.05, "Se rechaza H0", "No se rechaza H0"),
    ifelse(shapiro_sal$p.value < 0.05, "Se rechaza H0", "No se rechaza H0"),
    ifelse(bartlett_sal$p.value < 0.05, "Se rechaza H0", "No se rechaza H0"),
    ifelse(levene_sal$`Pr(>F)`[1] < 0.05, "Se rechaza H0", "No se rechaza H0"),
    ifelse(welch_sal$p.value < 0.05, "Se rechaza H0", "No se rechaza H0")
  )
)
colnames(resumen_c) <- c("Prueba", "Estadístico", "p-valor", "Decisión")

knitr::kable(resumen_c, row.names = FALSE, align = "c",
             caption = "Tabla 10. Todas las pruebas del punto c)")
Tabla 10. Todas las pruebas del punto c)
Prueba Estadístico p-valor Decisión
ANOVA (clásica) 29.8930 0.0000000 Se rechaza H0
Shapiro-Wilk (normalidad) 0.9732 0.3749000 No se rechaza H0
Bartlett (homocedasticidad) 20.0840 0.0000435 Se rechaza H0
Levene (homocedasticidad) 8.6750 0.0007020 Se rechaza H0
Welch (ANOVA robusto) 34.7270 0.0000001 Se rechaza H0

En conjunto, las cinco pruebas cuentan la misma historia: el pH sí genera diferencias significativas en la biomasa (ANOVA y Welch), los residuales del ANOVA clásico se comportan como normales (Shapiro), pero el supuesto de varianzas iguales no se cumple (Bartlett y Levene) — por eso la prueba de Welch, que no exige varianzas iguales, es la más confiable para la conclusión final. La prueba LSD muestra además que los tres niveles de pH (Bajo, Medio, Alto) difieren entre sí en biomasa.

Conclusión General

Los análisis desarrollados en este trabajo permitieron caracterizar la relación entre las propiedades del suelo y la producción de biomasa de la planta forrajera a partir de 45 muestras. El análisis exploratorio mostró diferencias importantes en la distribución y variabilidad de las variables. La biomasa presentó una variabilidad relativamente alta y cierta asimetría positiva, mientras que la salinidad mostró el comportamiento más homogéneo. El zinc y el potasio presentaron mayor variabilidad, con un grupo particular de valores bajos de zinc.

El análisis bivariado identificó al pH como la variable con mayor asociación con la biomasa, mostrando una relación positiva y un R² cercano a 0.86. El zinc también presentó una relación importante, aunque negativa, pero su asociación con el pH impide atribuirle un efecto independiente sobre la biomasa. Por otro lado, la salinidad y el potasio presentaron relaciones muy débiles y no mostraron capacidad relevante para explicar la variación de la biomasa en este conjunto de datos.

La categorización del pH mostró un aumento progresivo de la biomasa desde los niveles bajos hasta los altos, y el ANOVA confirmó que existen diferencias significativas entre las medias de los tres niveles de pH. Aunque los residuales cumplieron el supuesto de normalidad, las pruebas de Bartlett y Levene evidenciaron que las varianzas no son homogéneas, principalmente por la mayor variabilidad del grupo de pH alto. Por esta razón, la prueba de Welch resulta especialmente importante y confirmó que las diferencias entre los niveles de pH son significativas. Además, la prueba LSD mostró diferencias significativas entre los tres grupos.

En conjunto, los resultados muestran que, dentro de las variables estudiadas, el pH es la característica del suelo que presenta la asociación más fuerte con la producción de biomasa y que sus diferentes niveles se relacionan con diferencias significativas en la respuesta de la planta. Sin embargo, estos resultados representan asociaciones dentro de las muestras analizadas y no permiten afirmar que el pH, por sí solo, sea la causa directa del aumento de biomasa.

Punto 2 - Datos Moluscos

Contexto

La tasa de consumo de oxígeno es una medida crucial en la comprensión de los procesos fisiológicos de la gran mayoría de animales, lo cual incluye a un grupo de invertebrados que suele encontrarse más frecuentemente en ambientes marinos: los moluscos. Para poder estudiar la injerencia que puede llegar a tener la concentración/salinidad del agua marina y la morfofisiología de los moluscos en su tasa de consumo de oxígeno, se tomaron 48 datos de esta medida en 2 tipos distintos de molusco (A y B) y bajo 3 porcentajes distintos de concentración del agua (50, 75 y 100%).

1. Análisis exploratorio univariado

Previo a la realización de las pruebas estadísticas correspondientes, debe hacerse un proceso de observación y comprensión de cada uno de los indicadores estadísticos principales, tanto de tendencia central, como de dispersión y forma, principalmente. En este estudio se presentan 2 factores de tratamiento (concentración del agua y tipo de molusco) y una sola variable de respuesta. Por tanto, únicamente se realizará el análisis exploratorio de la respuesta, que es la cantidad de consumo de oxígeno, empleando las bibliotecas gt y htmltools para hacer visibles los cálculos de cada indicador:

resumen_m1 <- BD_moluscos %>%
  summarise(
    across(
      cons_o,
      list(
        n = ~sum(!is.na(.x)),
        Media = ~mean(.x, na.rm = TRUE),
        DE = ~sd(.x, na.rm = TRUE),
        Mediana = ~median(.x, na.rm = TRUE),
        Min = ~min(.x, na.rm = TRUE),
        Max = ~max(.x, na.rm = TRUE),
        CV = ~sd(.x, na.rm = TRUE) / mean(.x, na.rm = TRUE) * 100
      )
    )
  ) %>%
  pivot_longer(
    cols = everything(),
    names_to = c("Variable", "Estadistico"),
    names_pattern = "^(.*)_(n|Media|DE|Mediana|Min|Max|CV)$"
  ) %>%
  pivot_wider(
    names_from = Estadistico,
    values_from = value
  ) %>%
  select(Variable, n, Media, DE, Mediana, Min, Max, CV) %>%
  mutate(
    across(where(is.numeric), ~round(.x, 2))
    )
resumen_m1 <- resumen_m1 %>%
  rename(
    `Respuesta` = Variable,
    `n` = n,
    `Media` = Media,
    `DE` = DE,
    `Mediana` = Mediana,
    `Mínimo` = Min,
    `Máximo` = Max,
    `CV (%)` = CV
  )
resumen_m1 <- resumen_m1 %>%
  mutate(
      "Respuesta" = "Consumo oxígeno"
    )
tagList(
  tags$div(
    style = "font-size: 14px; color: #5f6368; margin-bottom: 10px;",
    "Tabla 1.Indicadores estadísticos principales de la medida de consumo de oxígeno"
  ),
  
  resumen_m1 %>%
    gt() %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()
  ) %>%
    cols_align(
      align = "center",
      columns = everything()
    )
)
Tabla 1.Indicadores estadísticos principales de la medida de consumo de oxígeno
Respuesta n Media DE Mediana Mínimo Máximo CV (%)
Consumo oxígeno 48 9.3 3.68 9.7 1.8 18.8 39.58

Ahora, además de los datos de los indicadores estadísticos de dispersión y tendencia central, se debe encontrar el nivel de simetría de los datos, para lo cual será necesario realizar un histograma de los datos, empleando la biblioteca ggplot2:

dist_o=ggplot(BD_moluscos, aes(x = cons_o)) +
  geom_histogram(bins = 10, fill = "#4682B4", color = "white") +
  labs(title = "Distribución de las medidas de consumo de oxígeno",
       x = "Consumo de oxígeno", y = "Frecuencia") +theme_minimal()+
  theme(axis.text.x = element_text(angle = 20, hjust = 1),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5))
ggplotly(dist_o)

Con esta gráfica, podría afirmarse que los datos son relativamente asímetricos de forma negativa, encontrando la mayor densidad de datos hacia 5,57. Sin embargo, al no ser una diferencia significativa, es posible decir, en cambio, que la distribución de los datos es simétrica.

2. Análisis exploratorio bivariado

Ahora, debido a que este estudio únicamente tiene en cuenta una variable de respuesta, es difícil realizar observaciones acerca de sus datos. Por tanto, se debe hacer un análisis bivariado, en el que primero se comparen uno a uno los tratamientos con la respuesta y, posteriormente, sea posible analizar la relación que existe entre los factores juntos y la respuesta, siguiendo las pautas del diseño factorial.

a. Relación de la concentración del agua con el consumo de oxígeno

En primera instancia, se realizará un análisis sobre el primer factor de tratamiento, correspondiente a la concentración del agua del mar. Para esto, se tabularán los datos de los indicadores de tendencia central de esta relación a continuación, empleando las bibliotecas Kable y KableExtra.

require(kableExtra)
tabla_cons_o_final <- BD_moluscos %>%
  summarise(
    `50 % - Media` = mean(cons_o[c_agua == 50], na.rm = TRUE),
    `50 % - DE` = sd(cons_o[c_agua == 50], na.rm = TRUE),
    `75 % - Media` = mean(cons_o[c_agua == 75], na.rm = TRUE),
    `75 % - DE` = sd(cons_o[c_agua == 75], na.rm = TRUE),
    `100 % - Media` = mean(cons_o[c_agua == 100], na.rm = TRUE),
    `100 % - DE` = sd(cons_o[c_agua == 100], na.rm = TRUE)
  ) %>%
  mutate(Respuesta = "Consumo de oxígeno") %>%
  select(Respuesta, `50 % - Media`, `50 % - DE`, `75 % - Media`, `75 % - DE`,
    `100 % - Media`, `100 % - DE`) %>%
  mutate(across(where(is.numeric), ~round(.x, 2)))
kable(tabla_cons_o_final,
  col.names = c("Respuesta", "Media", "DE", "Media", "DE", "Media", "DE"
  ), caption = "Tabla 2.Indicadores estadísticos principales entre concentración del agua y consumo de oxígeno",
  align = "ccccccc",
  escape = FALSE
) %>%
  add_header_above(
    c(
      " " = 1,
      "50 %" = 2,
      "75 %" = 2,
      "100 %" = 2
    ),
    align = "c"
  ) %>%
  add_header_above(
    c(" " = 1,
      "Concentración del agua" = 6), align = "c"
  ) %>%
  kable_styling(
    bootstrap_options = "bordered",
    full_width = FALSE,
    position = "center",
    font_size = 14
  ) %>%
  row_spec(0, bold = TRUE, align = "center") %>%
  column_spec(1, bold = TRUE, width = "3.5cm")
Tabla 2.Indicadores estadísticos principales entre concentración del agua y consumo de oxígeno
Concentración del agua
50 %
75 %
100 %
Respuesta Media DE Media DE Media DE
Consumo de oxígeno 12.25 3.2 6.99 2.8 8.67 3

Como adicional, se realiza un gráfico boxplot, utilizando ggplot2, para identificar las regiones de esta interacción que están más relacionadas, que se ven por las regiones que se traslapan, además de ver los datos atípicos que tenga la muestra.

dis_int_conc_cons <- ggplot(BD_moluscos, aes(x = as.factor(c_agua), y = cons_o)) +
  geom_boxplot(fill = c("#E76F51", "#2A9D8F", "#4361EE")) +
  labs(x = "Concentración del agua", y = "Consumo de oxígeno")

ggplotly(dis_int_conc_cons) %>%
layout(title = "",
       xaxis = list(title = "Concentración del agua",
                    showgrid = TRUE, gridcolor = "#D9D9D9"),
       yaxis = list(title = "Consumo de oxígeno",
                    showgrid = TRUE, gridcolor = "#D9D9D9"),
       paper_bgcolor = "#FFFFFF",
       plot_bgcolor = "#FFFFFF")

Estos resultados muestran que, en primer lugar, los datos obtenidos de las 3 concentraciones de agua estudiadas no difieren tanto entre sí, es decir, tienen desviaciones estándar muy similares. Además, puede afirmarse que , de los 3 tratamientos posibles, es en el 50% de concentración del agua de mar en donde los resultados varían más con respecto a las otras dos concentraciones estudiadas.

b. Relación del tipo de molusco con el consumo de oxígeno

Para este caso, se realizará el mismo procedimiento del inciso a, pero esta vez analizando la relación que existe entre el tipo de molusco y el consumo de oxígeno:

tabla_mol_final <- BD_moluscos %>%
  summarise(
    `A - Media` = mean(cons_o[molusco == "A"], na.rm = TRUE),
    `A - DE` = sd(cons_o[molusco == "A"], na.rm = TRUE),
    `B - Media` = mean(cons_o[molusco == "B"], na.rm = TRUE),
    `B - DE` = sd(cons_o[molusco == "B"], na.rm = TRUE),
  ) %>%
  mutate(Respuesta = "Consumo de oxígeno") %>%
  select(Respuesta, `A - Media`, `A - DE`, `B - Media`, `B - DE`) %>%
  mutate(across(where(is.numeric), ~round(.x, 2)))
kable(tabla_mol_final,
      col.names = c("Respuesta", "Media", "DE", "Media", "DE"), 
      caption = "Tabla 3.Indicadores estadísticos principales entre el tipo de molusco y el consumo de oxígeno",
      align = "ccccccc",
      escape = FALSE
) %>%
  add_header_above(
    c(
      " " = 1,
      "A" = 2,
      "B" = 2
    ),
    align = "c"
  ) %>%
  add_header_above(
    c(" " = 1,
      "Molusco" = 4), align = "c"
  ) %>%
  kable_styling(
    bootstrap_options = "bordered",
    full_width = FALSE,
    position = "center",
    font_size = 14
  ) %>%
  row_spec(0, bold = TRUE, align = "center") %>%
  column_spec(1, bold = TRUE, width = "3.5cm")
Tabla 3.Indicadores estadísticos principales entre el tipo de molusco y el consumo de oxígeno
Molusco
A
B
Respuesta Media DE Media DE
Consumo de oxígeno 10 3.27 8.61 4 
dis_int_mol_cons <- ggplot(BD_moluscos, aes(x = molusco, y = cons_o)) +
  geom_boxplot(fill = c("#E68F71", "#1368EE")) +
  labs(x = "Tipo de molusco", y = "Consumo de oxígeno")

ggplotly(dis_int_mol_cons) %>%
  layout(title = "",
         xaxis = list(title = "Tipo de molusco",
                      showgrid = TRUE, gridcolor = "#D9D9D9"),
         yaxis = list(title = "Consumo de oxígeno",
                      showgrid = TRUE, gridcolor = "#D9D9D9"),
         paper_bgcolor = "#FFFFFF",
         plot_bgcolor = "#FFFFFF")

En este caso, tanto el análisis de indicadores estadísticos como la representación mediante el boxplot indican que los datos de la interacción entre este factor y la respuesta, independientemente del tipo de molusco, indican una relación bastante estrecha entre sí, por lo que este factor no parece ser muy determinante en afectar la cantidad de consumo de oxígeno.

3. Análisis de varianza para la interacción entre variables (ANOVA)

Antes de realizar propiamente el análisis de varianzas, es necesario observar el gráfico que muestra la interacción entre todas las variables del estudio, que puede hacerse gracias a la biblioteca ggplot 2 y se muestra a continuación:

graf_int=ggplot(BD_moluscos, aes(x = c_agua, y = cons_o, colour = molusco))+geom_smooth()+geom_point()+labs(
  title = "Consumo de oxígeno de acuerdo con el tipo de molusco
  y la concentración del agua del mar",
  x = "Concentración del agua (%)",
  y = "Consumo de oxígeno",
  colour = "Tipo de molusco"
) +
  theme(plot.title = element_text(hjust = 0.5))
ggplotly(graf_int)

Ahora bien, una vez hecho todo este análisis exploratorio, se realizó un ANOVA para confirmar cuál de las variables del experimento está más relacionada con el consumo de oxígeno en los especímenes estudiados, obteniendo los siguientes resultados:

anova_mol=aov(cons_o~as.factor(c_agua)+molusco+as.factor(c_agua):molusco, data=BD_moluscos)
tabla=tidy(anova_mol)
tabla$term <- c("Concentración del agua", "Molusco", "Interacción", "Residuals")
tabla <- tabla %>%
  mutate(
    across(c(sumsq, meansq, statistic), ~ round(.x, 1)),
    p.value = ifelse(is.na(p.value), "", format.pval(p.value, digits = 3, eps = 0.001))
  )

kable(tabla,
      col.names = c("Fuente de variación", "Grados de libertad",
                    "Suma de cuadrados", "Cuadrado medio",
                    "Valor F", "p-valor"),
      caption = "Tabla 4. Anova de dos vías mostrando interacción entre  consumo de oxígeno, tipo de molusco y concentración del agua de mar",
      align = "c") %>%
  kable_styling(bootstrap_options = "bordered", full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE)
Tabla 4. Anova de dos vías mostrando interacción entre consumo de oxígeno, tipo de molusco y concentración del agua de mar
Fuente de variación Grados de libertad Suma de cuadrados Cuadrado medio Valor F p-valor
Concentración del agua 2 230.8 115.4 13.2 <0.001
Molusco 1 23.2 23.2 2.7 0.111
Interacción 2 15.4 7.7 0.9 0.424
Residuals 42 368.0 8.8 NA

Con la información de la Tabla 4 y la gráfica demostrando la interacción de los factores, pueden reafirmarse diversas ideas previas. Entre estas, se reafirma la poca diferencia existente entre los tratamientos hechos a partir del tipo de molusco, así como su poca incidencia en la cantidad de consumo de oxígeno. Asimismo, puede confirmarse que la variable que más incide en la cantidad de oxígeno consumido por el molusco es la concentración del agua de mar, siendo entonces que entre menor sea esta concentración, más elevada será la cantidad de oxígeno que consuma el animal.

Por último, para terminar de confirmar estas observaciones, se realizará un postANOVA, o prueba post-hoc, realizando la prueba de Shapiro-Wilk (con su respectivo histograma y gráfica Q-Q) y de homogeneidad de varianzas, las cuales se presentarán a continuación:

residuales_m <- residuals(anova_mol)
shapiro_residualesm <- shapiro.test(residuales_m)

tabla_shapiro_residm <- data.frame(
  Prueba = "Normalidad de residuales (Anova interacción concentración del agua y tipo de molusco)",
  Estadistico_W = round(unname(shapiro_residualesm$statistic), 3),
  Valor_p = round(shapiro_residualesm$p.value, 4),
  Normalidad = ifelse(shapiro_residualesm$p.value > 0.05, "No se rechaza", "Se rechaza")
)

kable(tabla_shapiro_residm,
      col.names = c("Prueba", "Estadístico W", "Valor p", "Normalidad (α = 0.05)"),
      caption = "Prueba de Shapiro-Wilk sobre los residuales del modelo",
      row.names = FALSE,
      align = "lccc")
Prueba de Shapiro-Wilk sobre los residuales del modelo
Prueba Estadístico W Valor p Normalidad (α = 0.05)
Normalidad de residuales (Anova interacción concentración del agua y tipo de molusco) 0.958 0.0857 No se rechaza 
par(mfrow = c(1,2))
hist(residuales_m, main = "Histograma de residuales", xlab = "Residuales", col = "lightblue")
qqnorm(residuales_m); qqline(residuales_m, col = "red")

levene_conc <- leveneTest(cons_o ~ as.factor(c_agua), data = BD_moluscos)

tabla_levenem <- as.data.frame(levene_conc)
rownames(tabla_levenem) <- NULL
tabla_levenem$`Pr(>F)` <- format.pval(tabla_levenem$`Pr(>F)`, digits = 3, eps = 0.001)
tabla_levenem$`F value` <- round(tabla_levenem$`F value`, 3)

kable(tabla_levenem,
      col.names = c("GL", "F", "Valor p"),
      caption = "Prueba de homogeneidad de varianzas de concentración del agua entre hábitats",
      align = "ccc",
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = TRUE,
                position = "center",
                font_size = 14)
Prueba de homogeneidad de varianzas de concentración del agua entre hábitats
GL F Valor p
2 0.217 0.806
45 NA NA

Estas pruebas demuestran que no existe variabilidad grande entre los grupos de los factores, ni siquiera en los tratamientos de mayor diversidad como lo son los de la concentración del agua del mar. Sin embargo, esto sólo dice que las varianzas entre estos grupos son muy similares, no habla como tal de la injerencia que tengan sobre la variable de respuesta.

Para poder conocer la respuesta a la última pregunta, es necesario realizar un test LSD y graficar estos resultados, comparando nuevamente a los valores de consumo de oxígeno con los de concentración del agua de mar.

require(agricolae)
lsd_mol <- LSD.test(anova_mol, "as.factor(c_agua)", p.adj = "none", console = FALSE)
plot(lsd_mol,
     main = "Comparación de consumo de oxígeno por concentración del agua (LSD)",
     xlab = "Concentración del agua", ylab = "Consumo de oxígeno")

En este gráfico se observa entonces que sí existe una variación entre los tratamientos del factor de concentración del agua de mar, donde los tratamientos del 75 y 100% de concentración tienen el mismo grupo asignado (b) al tener medias similares; mientras que, el tratamiento de 50% de concentración de agua de mar queda aislado en el grupo a por tener una media mucho mayor a la del resto de tratamientos, tal como el análisis exploratorio y el ANOVA habían previsto.

Conclusiones

Tras las pruebas realizadas y el análisis de todos los datos estudiados junto a sus interacciones, es posible concluir que:

  1. No existe una relación real entre el tipo de molusco y un cambio en la cantidad de consumo de oxígeno del animal.
  2. Tampoco existen cambios drásticos al estudiar en conjunto, como interacción, el tipo de molusco y la concentración del agua de mar con la cantidad de consumo de oxígeno.
  3. Existe una fuerte relación entre una baja concentración del agua de mar y una alta cantidad de consumo de oxígeno, haciendo de este el único factor incidente para la variable de respuesta estudiada. Se esperaría entonces que los moluscos que viven en zonas con concentración de agua marina menores al 50% tiendan a necesitar una mucho mayor cantidad de oxígeno en comparación a otros moluscos que vivan en zonas con concentraciones mayores al 50%, aproximadamente.

Punto 3 - Datos Biodiversidad

Contexto

La pérdida y transformación de hábitats naturales es una de las principales causas de declive de la biodiversidad a nivel global, y los anfibios son considerados uno de los grupos más sensibles a estos cambios, debido a su piel permeable, sus ciclos de vida dependientes de humedad y su baja capacidad de dispersión frente a paisajes fragmentados. En este estudio se evaluó el efecto del uso del suelo sobre la biodiversidad de anfibios, comparando cuatro tipos de hábitat que representan un gradiente de intervención antrópica: bosque primario, bosque secundario, sistema silvopastoril y potrero.

Se establecieron 52 parcelas de muestreo distribuidas equitativamente entre los cuatro hábitats (13 parcelas por hábitat). En cada parcela se registraron tres variables: la riqueza de especies, el índice de diversidad de Shannon-Wiener y la altitud del sitio de muestreo.

A. Analisis Exploratorio Univariado

Antes de empezar con las pruebas estadisticas, es necesario conocer el comportamiento de cada variable registrada en el estudio. Esta primera evaluación, llamada análisis exploratorio univariado, aún no busca responder las preguntas centrales del estudio, sino ver como se distribuyen los datos, identificar valores atípicos o patrones inusuales, y verificar que las variables cumplen ciertos los supuestos, verificando que las conclusiones posteriores estén respaldadas por una comprensión real de la información.

BD_biodiversidad$Habitat <- factor(
  BD_biodiversidad$Habitat,
  levels = c("Bosque primario", "Bosque secundario",
             "Sistema silvopastoril", "Potrero")
)

Para variables registradas: riqueza de especies, índice de diversidad de Shannon-Wiener y altitud, se calcularon las medidas de resumen típicas: la media, la mediana, la desviación estándar y el coeficiente de variación.

La riqueza de especies presentó un promedio de 9.73 especies por parcela, con una desviación estándar de 4.37 y un rango de 2 hasta 21 especies observadas; su coeficiente de variación, de 44.9%, fue el más alto de las tres variables.

El índice de Shannon tuvo un promedio de 1.76 y una desviación estándar de 0.61, con un coeficiente de variación de 34.4%; su media y su mediana (1.85) resultaron muy cercanas entre sí, mostrando una distribución sin asimetrias marcadas.

Finalmente, la altitud tuvo un promedio de 1243.9 m.s.n.m., con una desviación estándar de 172.7 metros y, en comparación a las otras variables, un coeficiente de variación mucho más bajo, de apenas 13.9%, lo que sugiere que esta variable es mucho más homogénea.

library(tidyr)

resumen_univar <- BD_biodiversidad %>%
  summarise(across(c(Riqueza, Shannon, Altitud),
                   list(Media = ~mean(.x), DE = ~sd(.x),
                        Mediana = ~median(.x),
                        Min = ~min(.x), Max = ~max(.x),
                        CV = ~sd(.x)/mean(.x)*100))) %>%
  pivot_longer(cols = everything(),
               names_to = c("Variable", "Estadistico"),
               names_sep = "_") %>%
  pivot_wider(names_from = Estadistico, values_from = value) %>%
  mutate(across(where(is.numeric), ~round(.x, 2)))

resumen_univar

Estas distribuciones se examinaron gráficamente mediante histogramas realizados con la librería ggplot2.

El histograma de la riqueza mostró una forma simétrica, con la mayor concentración de parcelas entre 3 y 15 especies, coherente con los valores de media y mediana obtenidos previamente.

El histograma del índice de Shannon evidenció una alta concentración de observaciones alrededor de un valor cercano a 2, disminuyendo de forma gradual hacia ambos extremos, lo que refuerza la idea de una distribución sin asimetrías importantes.

El comportamiento de la altitud fue distinto: mostró una distribución con poca frecuencia en la zona central y una acumulación en los extremos, es decir, una forma que se aparta de la campana de distribución normal. Esto no es un error en los datos, sino que está determinada por la ubicación de cada tipo de hábitat, de modo que el histograma conjunto de las 52 parcelas está mezclando, en realidad, cuatro subpoblaciones distintas con rangos altitudinales diferentes.

ggplot(BD_biodiversidad, aes(x = Riqueza)) +
  geom_histogram(bins = 10, fill = "#2E8B57", color = "white") +
  labs(title = "Distribución de la riqueza de especies",
       x = "Riqueza (N° especies)", y = "Frecuencia") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5))

ggplot(BD_biodiversidad, aes(x = Shannon)) +
  geom_histogram(bins = 10, fill = "#4682B4", color = "white") +
  labs(title = "Distribucion del indice de Shannon-Wiener",
       x = "Indice de Shannon", y = "Frecuencia") +theme_minimal()+
  theme(axis.text.x = element_text(angle = 20, hjust = 1),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5))

ggplot(BD_biodiversidad, aes(x = Altitud)) +
  geom_histogram(bins = 10, fill = "#B8860B", color = "white") +
  labs(title = "Distribucion de la altitud",
       x = "Altitud (m.s.n.m.)", y = "Frecuencia") +  theme_minimal()+
  theme(axis.text.x = element_text(angle = 20, hjust = 1),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5))

Para respaldar estas observaciones, se aplicó la prueba de Shapiro-Wilk a cada una de las tres variables, cuya hipótesis nula establece que los datos provienen de una distribución normal.

Para la riqueza de especies se obtuvo un valor p de 0.053, en el límite del umbral convencional de significancia de 0.05, por lo que no se rechaza la normalidad, aunque el margen es estrecho.

Para el índice de Shannon se obtuvo un valor p de 0.389, claramente compatible con una distribución normal

Para la altitud se obtuvo un valor p de 0.025, que lleva a rechazar el supuesto de normalidad. Esto confirma lo que ya habiamos mencionado: la altitud, al estar asociada estrechamente al tipo de hábitat, no puede tratarse como una variable homogénea.

library(knitr)

shapiro_riqueza <- shapiro.test(BD_biodiversidad$Riqueza)
shapiro_shannon <- shapiro.test(BD_biodiversidad$Shannon)
shapiro_altitud <- shapiro.test(BD_biodiversidad$Altitud)

tabla_shapiro <- data.frame(
  Variable = c("Riqueza", "Shannon", "Altitud"),
  Estadistico_W = c(shapiro_riqueza$statistic,
                    shapiro_shannon$statistic,
                    shapiro_altitud$statistic),
  Valor_p = c(shapiro_riqueza$p.value,
              shapiro_shannon$p.value,
              shapiro_altitud$p.value)
)

tabla_shapiro$Estadistico_W <- round(tabla_shapiro$Estadistico_W, 3)
tabla_shapiro$Valor_p <- round(tabla_shapiro$Valor_p, 4)

kable(tabla_shapiro,
      col.names = c("Variable", "Estadístico W", "Valor p"),
      caption = "Prueba de Shapiro-Wilk para las variables originales")
Prueba de Shapiro-Wilk para las variables originales
Variable Estadístico W Valor p
Riqueza 0.956 0.0532
Shannon 0.976 0.3885
Altitud 0.949 0.0252

B. Analisis Exploratorio Bivariado

El siguiente paso consiste en explorar cómo se relacionan entre sí, y en particular cómo se comportan la riqueza y el índice de Shannon cuando se separan según el tipo de hábitat. Este análisis bivariado permite anticipar visualmente si existe una asociación entre el hábitat y la biodiversidad, antes de someter esa pregunta a una prueba formal de hipótesis; y tambien permite aclarar observaciones del análisis univariado, como el comportamiento de la altitud.

Para comparar la riqueza de especies entre los cuatro tipos de hábitat se construyó un boxplot, complementado con los puntos individuales de cada parcela en el diagrama, de manera que fuera posible observar tanto el resumen estadístico de cada grupo como la dispersión real de los datos originales. El gráfico mostró un patrón descendente claro: las cajas correspondientes al bosque primario se ubicaron en los valores más altos de riqueza, seguidas por el bosque secundario, luego el sistema silvopastoril, y finalmente el potrero, con las cajas más bajas de las cuatro.

El mismo ejercicio se repitió para el índice de Shannon, obteniéndose un patrón prácticamente idéntico: las cajas se ordenaron de la misma manera y con una separación similar entre grupos, lo que indica que ambas medidas de biodiversidad tienen un mismo patrón con respecto al habitat.

ggplot(BD_biodiversidad, aes(x = Habitat, y = Riqueza, fill = Habitat)) +
  geom_boxplot(show.legend = FALSE) +
  geom_jitter(width = 0.15, alpha = 0.4) +
  labs(title = "Riqueza de especies segun tipo de habitat",
       x = "Habitat", y = "Riqueza (# especies)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1),
        legend.position = "none",
plot.title = element_text(hjust = 0.5))

ggplot(BD_biodiversidad, aes(x = Habitat, y = Shannon, fill = Habitat)) +
  geom_boxplot(show.legend = FALSE) +
  geom_jitter(width = 0.15, alpha = 0.4) +
  labs(title = "Indice de Shannon-Wiener segun tipo de habitat",
       x = "Habitat", y = "Indice de Shannon") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1),
        legend.position = "none",
plot.title = element_text(hjust = 0.5))

Para reespaldar estas observaciones se calcularon las medias y desviaciones estándar de riqueza, Shannon y altitud para cada uno de los cuatro hábitats. El bosque primario presentó una riqueza promedio de 14.8 especies y un índice de Shannon promedio de 2.47; el bosque secundario, 11.1 especies y un Shannon de 1.97; el sistema silvopastoril, 8.1 especies y un Shannon de 1.60; y el potrero, apenas 5.0 especies y un Shannon de 1.00.

Es decir, el bosque primario tiene en promedio casi tres veces más especies que el potrero, y un índice de diversidad más de dos veces mayor. Cabe destacar que las desviaciones estándar del índice de Shannon resultaron similares para cada habitat (entre 0.22 y 0.38).

library(kableExtra)

tabla_habitat <- BD_biodiversidad %>%
  group_by(Habitat) %>%
  summarise(Media_Riqueza = mean(Riqueza), DE_Riqueza = sd(Riqueza),
            Media_Shannon = mean(Shannon), DE_Shannon = sd(Shannon),
            Media_Altitud = mean(Altitud), DE_Altitud = sd(Altitud)) %>%
  mutate(across(where(is.numeric), ~round(.x, 2)))

kable(tabla_habitat,
      col.names = c("Hábitat", "Media", "DE", "Media", "DE", "Media", "DE"),
      caption = "Estadísticos descriptivos de riqueza, Shannon y altitud por tipo de hábitat",
      align = "lcccccc") %>%
  add_header_above(c(" " = 1, "Riqueza" = 2, "Shannon" = 2, "Altitud" = 2)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE,
                position = "center",
                font_size = 14)
Estadísticos descriptivos de riqueza, Shannon y altitud por tipo de hábitat
Riqueza
Shannon
Altitud
Hábitat Media DE Media DE Media DE
Bosque primario 14.77 2.74 2.47 0.38 1453.19 77.50
Bosque secundario 11.08 3.01 1.97 0.22 1326.15 40.60
Sistema silvopastoril 8.08 2.50 1.60 0.28 1151.48 97.89
Potrero 5.00 1.29 1.00 0.24 1044.75 50.32

Se exploró la relación entre la altitud y el índice de Shannon mediante un gráfico de dispersión, coloreando cada punto según su hábitat, junto con una línea de tendencia sobre el conjunto de datos.

Este gráfico permite explicar el histograma de la altitud que vimos en el análisis univariado: los puntos no se distribuyeron al azar, sino que se agruparon en franjas diferenciadas según cada hábitat. El bosque primario ocupando las altitudes más elevadas y el potrero las más bajas.

El cálculo de la correlación entre ambas variables confirmó una asociación positiva fuerte y altamente significativa (r = 0.75, p < 0.001). Esto es una correlación, y no una causalidad, puesto que en realidad la altitud está confundida con el tipo de hábitat: dado que cada tipo de uso del suelo se concentra en una altitud particular dentro de la reserva, no es posible separar si la diversidad cambia por efecto de la altitud o por efecto del hábitat.

ggplot(BD_biodiversidad, aes(x = Altitud, y = Shannon, color = Habitat)) +
  geom_point(size = 2) +
  geom_smooth(aes(group = 1), method = "lm", se = TRUE, color = "black", linetype = "dashed") +
  labs(title = "Relacion entre altitud y diversidad de Shannon",
       x = "Altitud (m.s.n.m.)", y = "Indice de Shannon") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1),
plot.title = element_text(hjust = 0.5))

correlacion <- cor.test(BD_biodiversidad$Altitud, BD_biodiversidad$Shannon, method = "pearson")

tabla_correlacion <- data.frame(
  Variables = "Altitud vs. Shannon",
  r = round(unname(correlacion$estimate), 3),
  IC_95_inf = round(unname(correlacion$conf.int[1]), 3),
  IC_95_sup = round(unname(correlacion$conf.int[2]), 3),
  Valor_p = format.pval(correlacion$p.value, digits = 3, eps = 0.001)
)

kable(tabla_correlacion,
      col.names = c("Variables", "r de Pearson", "IC 95% inferior", "IC 95% superior", "Valor p"),
      caption = "Correlación entre altitud y diversidad de Shannon",
      row.names = FALSE,
      align = "lcccc") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE,
                position = "center",
                font_size = 14)
Correlación entre altitud y diversidad de Shannon
Variables r de Pearson IC 95% inferior IC 95% superior Valor p
Altitud vs. Shannon 0.752 0.603 0.851 <0.001

C. Analisis de varianza (ANOVA) y comparaciones post-hoc

Los dos análisis anteriores sugirieron que el tipo de hábitat está asociado a diferencias en la diversidad de anfibios, sin embargo, necesitamos responder esa pregunta de manera formal. Para determinar si las diferencias entre hábitats son suficientemente grandes como para no atribuirse simplemente al azar, se ajustó un análisis de varianza (ANOVA) de una vía, en el cual la variable respuesta fue el índice de diversidad de Shannon y el único factor fue el tipo de hábitat, con sus cuatro niveles: bosque primario, bosque secundario, sistema silvopastoril y potrero.

En este modelo se descompone la variabilidad total observada en el índice de Shannon en dos componentes: la variabilidad entre las medias de los cuatro hábitats, y la variabilidad dentro de cada hábitat. Si la primera resulta considerablemente mayor que la segunda, se interpreta como evidencia de que el hábitat sí genera diferencias reales en la diversidad, más allá de la variación natural que cabría esperar entre parcelas similares.

El resultado del ANOVA mostró que la suma de cuadrados asociada al hábitat fue de 14.86, frente a una suma de cuadrados residual de apenas 3.92, es decir, la variabilidad explicada por las diferencias entre hábitats resultó casi cuatro veces mayor que la variabilidad no explicada dentro de los grupos. Esta relación se resume en el estadístico F, que en este caso fue de 60.61. El valor p asociado a este estadístico fue extremadamente pequeño (p < 0.001), muy por debajo del umbral convencional de 0.05, por lo que se concluye que el tipo de hábitat tiene un efecto altamente significativo sobre el índice de diversidad de Shannon. En otras palabras, existe evidencia estadística sólida de que al menos uno de los cuatro hábitats presenta una diversidad promedio distinta a la de los demás, resultado que respalda de manera formal el patrón que ya se había insinuado en los boxplots del literal anterior.

modelo_shannon <- aov(Shannon ~ Habitat, data = BD_biodiversidad)
tabla_anova <- as.data.frame(summary(modelo_shannon)[[1]])

tabla_anova$Df <- as.integer(tabla_anova$Df)
tabla_anova$`Sum Sq` <- round(tabla_anova$`Sum Sq`, 3)
tabla_anova$`Mean Sq` <- round(tabla_anova$`Mean Sq`, 3)
tabla_anova$`F value` <- round(tabla_anova$`F value`, 2)
tabla_anova$`Pr(>F)` <- format.pval(tabla_anova$`Pr(>F)`, digits = 3, eps = 0.001)

kable(tabla_anova,
      col.names = c("GL", "Suma de cuadrados", "Cuadrado medio", "F", "Valor p"),
      caption = "Tabla ANOVA: Shannon en función del hábitat",
      align = "ccccc") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE,
                position = "center",
                font_size = 14)
Tabla ANOVA: Shannon en función del hábitat
GL Suma de cuadrados Cuadrado medio F Valor p
Habitat 3 14.862 4.954 60.61 <0.001
Residuals 48 3.923 0.082 NA NA

Antes de dar por válida esta conclusión es necesario comprobar que el modelo cumple los dos supuestos sobre los cuales se sustenta la validez del ANOVA: la normalidad de los residuales y la homogeneidad de varianzas entre los grupos.

Estos supuestos se evalúan sobre la diferencia entre cada valor observado y la media de su propio grupo de hábitat, ya que es la parte no explicada por el factor la que debe comportarse como ruido aleatorio para que las pruebas de significancia del ANOVA sean confiables.

La normalidad de los residuales se evaluó mediante la prueba de Shapiro-Wilk, complementada con un histograma y un gráfico de cuantiles teóricos (QQ-plot). El resultado de la prueba (W = 0.978, p = 0.4625) no permite rechazar la normalidad, y el gráfico QQ mostró que los puntos se ajustan de forma cercana a la línea teórica, sin desviaciones sistemáticas en los extremos.

La homogeneidad de varianzas se evaluó mediante la prueba de Levene centrada en la mediana, la cual tampoco encontró evidencia de que las varianzas difirieran entre los cuatro hábitats (F = 1.277, p = 0.293). Al cumplirse ambos supuestos, se confirma que el ANOVA es válido y que sus conclusiones son confiables.

residuales <- residuals(modelo_shannon)
shapiro_residuales <- shapiro.test(residuales)

tabla_shapiro_resid <- data.frame(
  Prueba = "Normalidad de residuales (ANOVA Shannon ~ Hábitat)",
  Estadistico_W = round(unname(shapiro_residuales$statistic), 3),
  Valor_p = round(shapiro_residuales$p.value, 4),
  Normalidad = ifelse(shapiro_residuales$p.value > 0.05, "No se rechaza", "Se rechaza")
)

kable(tabla_shapiro_resid,
      col.names = c("Prueba", "Estadístico W", "Valor p", "Normalidad (α = 0.05)"),
      caption = "Prueba de Shapiro-Wilk sobre los residuales del modelo",
      row.names = FALSE,
      align = "lccc")
Prueba de Shapiro-Wilk sobre los residuales del modelo
Prueba Estadístico W Valor p Normalidad (α = 0.05)
Normalidad de residuales (ANOVA Shannon ~ Hábitat) 0.978 0.4625 No se rechaza 
par(mfrow = c(1,2))
hist(residuales, main = "Histograma de residuales", xlab = "Residuales", col = "lightblue")
qqnorm(residuales); qqline(residuales, col = "red")

par(mfrow = c(1,1))
levene_shannon <- leveneTest(Shannon ~ Habitat, data = BD_biodiversidad)

tabla_levene <- as.data.frame(levene_shannon)
rownames(tabla_levene) <- NULL
tabla_levene$`Pr(>F)` <- format.pval(tabla_levene$`Pr(>F)`, digits = 3, eps = 0.001)
tabla_levene$`F value` <- round(tabla_levene$`F value`, 3)

kable(tabla_levene,
      col.names = c("GL", "F", "Valor p"),
      caption = "Prueba de Levene: homogeneidad de varianzas de Shannon entre hábitats",
      align = "ccc",
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = TRUE,
                position = "center",
                font_size = 14)
Prueba de Levene: homogeneidad de varianzas de Shannon entre hábitats
GL F Valor p
3 1.277 0.293
48 NA NA

Dado que el ANOVA fue significativo, lo siguiente es identificar cuáles de los cuatro hábitats difieren entre sí, ya que la prueba global únicamente indica que existe al menos una diferencia, sin señalar entre qué grupos ocurre.

Para ello se aplicó la prueba de comparaciones múltiples de la diferencia mínima significativa (LSD), la cual compara las medias de cada par posible de hábitats utilizando la variabilidad residual estimada del ANOVA. Se optó por esta prueba, sin ajuste adicional por comparaciones múltiples, dado que el número de grupos a comparar es reducido (cuatro hábitats, seis comparaciones posibles) y el ANOVA global ya había arrojado un resultado significativo.

El resultado de esta prueba asignó una letra distinta a cada uno de los cuatro hábitats, lo que significa que ningún par de hábitats comparte una misma letra y, por lo tanto, que las seis comparaciones posibles resultaron estadísticamente significativas. Es decir, no solo el bosque primario difiere del potrero, como podría esperarse al comparar los dos extremos del gradiente, sino que cada hábitat difiere significativamente de todos los demás, incluyendo pares consecutivos como bosque secundario frente a sistema silvopastoril, o sistema silvopastoril frente a potrero.

lsd_shannon <- LSD.test(modelo_shannon, "Habitat", p.adj = "none", console = FALSE)
plot(lsd_shannon,
     main = "Comparación de medias de Shannon por hábitat (LSD)",
     xlab = "Hábitat", ylab = "Índice de Shannon")

Conclusión General

Los analisis desarrollados en este trabajo permiten responder consistentemente a la pregunta sobre el efecto del uso del suelo en la biodiversidad de anfibios de la reserva forestal. El análisis exploratorio univariado mostró que la riqueza de especies y el índice de Shannon se comportan de manera parecida a una distribución normal, mientras que la altitud presentó un comportamiento distinto, explicado por su estrecha asociación con el tipo de hábitat más que por un problema en los datos mismos. El análisis bivariado permitió observar que tanto la riqueza como la diversidad de Shannon disminuyen de manera progresiva a lo largo del gradiente de intervención antrópica, y que esta misma asociación con el hábitat es la que explica el comportamiento particular de la altitud detectado previamente, al estar esta variable confundida con el tipo de uso del suelo y no actuar como un factor independiente.

El análisis de varianza confirmó de manera formal lo que los gráficos nos mostraban: el tipo de hábitat tiene un efecto altamente significativo sobre el índice de diversidad de Shannon, y este resultado es plenamente confiable dado que el modelo cumplió satisfactoriamente los supuestos de normalidad de los residuales y de homogeneidad de varianzas. Las comparaciones post-hoc de LSD, a su vez, revelaron que esta diferencia no se limita a una separación entre hábitats extremos, sino que constituye un gradiente completo en el que cada tipo de hábitat presenta una diversidad significativamente distinta a la de todos los demás, incluyendo pares de hábitats consecutivos y aparentemente similares en su estructura.

En conjunto, estos resultados aportan evidencia sólida de que la transformación del bosque primario, incluso cuando conserva cierta cobertura arbórea como ocurre en el bosque secundario o en el sistema silvopastoril, tiene un costo medible y progresivo sobre la comunidad de anfibios de la reserva. Esto tiene implicaciones relevantes para el manejo y la conservación del área de estudio: la protección del bosque primario resulta insustituible para mantener los niveles más altos de biodiversidad, pero al mismo tiempo, los resultados sugieren que el mantenimiento de sistemas con algún grado de cobertura arbórea, como el silvopastoril, cumple un papel intermedio valioso frente a la alternativa de un potrero completamente desprovisto de árboles, y podría constituir una estrategia razonable en zonas donde la conservación estricta del bosque primario no sea viable.