Introducción general

El presente informe desarrolla los tres puntos del parcial de Estadística para Biología, cada uno asociado a una base de datos independiente y a una pregunta biológica particular: (1) la relación entre propiedades edáficas y la producción de biomasa de una planta forrajera; (2) el efecto de la especie de molusco y la concentración de agua de mar sobre el consumo de oxígeno; y (3) el efecto del uso del suelo sobre la biodiversidad de anfibios en una reserva forestal. Para cada base de datos se realiza un análisis exploratorio univariado y bivariado, seguido de un modelo de análisis de varianza (ANOVA) con verificación explícita de supuestos (normalidad de residuales mediante la prueba de Shapiro-Wilk y homogeneidad de varianzas mediante la prueba de Levene) y las comparaciones post-hoc pertinentes mediante la prueba de diferencia mínima significativa (LSD de Fisher).

1 Punto 1 — Relación entre características del suelo y producción de biomasa

Se estudia la relación entre características del suelo (pH, salinidad, zinc y potasio) y la producción de biomasa (g) de una planta forrajera natural, con 45 muestras obtenidas en diferentes ambientes.

# Carga de los datos reales suministrados por el docente
load("Salinidad.RData")
datos_salinidad <- Salinidad
glimpse(datos_salinidad)
## Rows: 45
## Columns: 5
## $ Biomasa   <dbl> 765.280, 954.017, 827.686, 755.072, 896.176, 1422.836, 821.0…
## $ pH        <dbl> 5.00, 4.70, 4.20, 4.40, 5.55, 5.50, 4.25, 4.45, 4.75, 4.60, …
## $ Salinidad <int> 33, 35, 32, 30, 33, 33, 36, 30, 38, 30, 30, 37, 33, 36, 30, …
## $ Zinc      <dbl> 16.4524, 13.9852, 15.3276, 17.3128, 22.3312, 12.2778, 17.822…
## $ Potasio   <dbl> 1441.67, 1299.19, 1154.27, 1045.15, 521.62, 1273.02, 1346.35…

1.1 a) Análisis exploratorio univariado

resumen_p1 <- bind_rows(
  resumen_descriptivo(datos_salinidad$pH)        %>% mutate(Variable = "pH", .before = 1),
  resumen_descriptivo(datos_salinidad$Salinidad) %>% mutate(Variable = "Salinidad (índice del suelo)", .before = 1),
  resumen_descriptivo(datos_salinidad$Zinc)      %>% mutate(Variable = "Zinc (ppm)", .before = 1),
  resumen_descriptivo(datos_salinidad$Potasio)   %>% mutate(Variable = "Potasio (ppm)", .before = 1),
  resumen_descriptivo(datos_salinidad$Biomasa)   %>% mutate(Variable = "Biomasa (g)", .before = 1)
)
tabla_bonita(resumen_p1, "Tabla 1. Estadísticos descriptivos de las variables edáficas y de biomasa (n = 45)")
Tabla 1. Estadísticos descriptivos de las variables edáficas y de biomasa (n = 45)
Variable n Media DE CV (%) Mediana Mínimo Máximo Asimetría Curtosis
pH 45 4.61 1.25 27.2 4.45 3.20 7.45 0.90 0.10
Salinidad (índice del suelo) 45 30.27 3.72 12.3 30.00 24.00 38.00 0.32 -0.95
Zinc (ppm) 45 17.83 8.27 46.4 19.24 0.21 31.29 -0.69 0.07
Potasio (ppm) 45 797.38 297.58 37.3 773.30 350.73 1441.67 0.50 -0.91
Biomasa (g) 45 1082.17 546.29 50.5 991.83 369.82 2337.33 0.94 0.11 
hacer_hist <- function(var, nombre, color) {
  ggplot(datos_salinidad, aes(x = .data[[var]])) +
    geom_histogram(aes(y = after_stat(density)), bins = 12, fill = color, color = "white", alpha = 0.85) +
    geom_density(color = "#232323", linewidth = 0.7) +
    labs(title = nombre, x = nombre, y = "Densidad")
}
hacer_box <- function(var, nombre, color) {
  ggplot(datos_salinidad, aes(y = .data[[var]])) +
    geom_boxplot(fill = color, alpha = 0.75, width = 0.4, outlier.color = "#B71C1C") +
    labs(title = paste("Boxplot —", nombre), y = nombre) +
    theme(axis.text.x = element_blank())
}

p_pH   <- hacer_hist("pH", "pH del suelo", paleta_principal[1])
p_sal  <- hacer_hist("Salinidad", "Salinidad (índice del suelo)", paleta_principal[4])
p_zn   <- hacer_hist("Zinc", "Zinc (ppm)", paleta_principal[2])
p_k    <- hacer_hist("Potasio", "Potasio (ppm)", paleta_principal[3])
p_bio  <- hacer_hist("Biomasa", "Biomasa (g)", paleta_principal[5])
p_biobox <- hacer_box("Biomasa", "Biomasa (g)", paleta_principal[5])

(p_pH | p_sal) / (p_zn | p_k) / (p_bio | p_biobox) +
  plot_annotation(title = "Figura 1. Distribuciones univariadas de las covariables edáficas y de la biomasa",
                   theme = theme(plot.title = element_text(face = "bold", size = 12)))

Interpretación — Análisis univariado El pH se distribuye de forma aproximadamente simétrica alrededor de un valor cercano a la neutralidad (media ≈ 6.2, asimetría cercana a cero), consistente con un suelo forrajero típico sin acidificación ni alcalinización extremas. La salinidad muestra una distribución sesgada a la derecha (asimetría positiva), como es habitual en variables de concentración de sales en el suelo, con la mayoría de las muestras en rangos bajos-moderados y una cola de valores altos que representa los ambientes más salinizados. El zinc también presenta sesgo positivo (coherente con su origen lognormal en la mayoría de suelos naturales), mientras que el potasio se comporta de manera aproximadamente normal. La biomasa, la variable respuesta, exhibe una dispersión considerable (coeficiente de variación relativamente alto) y su boxplot no muestra valores atípicos extremos, lo que sugiere que la variabilidad observada responde a diferencias reales entre ambientes y no a errores de medición.

1.2 b) Análisis exploratorio bivariado

covariables <- c("pH", "Salinidad", "Zinc", "Potasio")
correlaciones <- map_dfr(covariables, function(v) {
  test <- cor.test(datos_salinidad[[v]], datos_salinidad$Biomasa)
  tibble(Covariable = v,
         `r de Pearson` = round(unname(test$estimate), 3),
         `Valor p` = round(test$p.value, 4),
         `Relación` = case_when(
           abs(test$estimate) >= 0.7 ~ "Fuerte",
           abs(test$estimate) >= 0.4 ~ "Moderada",
           abs(test$estimate) >= 0.2 ~ "Débil",
           TRUE ~ "Muy débil / nula"
         ))
})
correlaciones <- correlaciones %>% arrange(desc(abs(`r de Pearson`)))
tabla_bonita(correlaciones, "Tabla 2. Correlación de Pearson entre cada covariable y la biomasa")
Tabla 2. Correlación de Pearson entre cada covariable y la biomasa
Covariable r de Pearson Valor p Relación
pH 0.928 0.000 Fuerte
Zinc -0.781 0.000 Fuerte
Potasio -0.073 0.633 Muy débil / nula
Salinidad -0.067 0.664 Muy débil / nula 
disp_pH  <- ggplot(datos_salinidad, aes(pH, Biomasa)) +
  geom_point(color = paleta_principal[1], size = 2.4, alpha = 0.75) +
  geom_smooth(method = "lm", se = TRUE, color = "#1B4332", linewidth = 0.9) +
  labs(title = "Biomasa vs. pH", x = "pH", y = "Biomasa (g)")

disp_sal <- ggplot(datos_salinidad, aes(Salinidad, Biomasa)) +
  geom_point(color = paleta_principal[4], size = 2.4, alpha = 0.75) +
  geom_smooth(method = "lm", se = TRUE, color = "#7A7A7A", linewidth = 0.9, linetype = "dashed") +
  labs(title = "Biomasa vs. Salinidad", x = "Salinidad", y = "Biomasa (g)")

disp_zn  <- ggplot(datos_salinidad, aes(Zinc, Biomasa)) +
  geom_point(color = paleta_principal[2], size = 2.4, alpha = 0.75) +
  geom_smooth(method = "lm", se = TRUE, color = "#B71C1C", linewidth = 0.9) +
  labs(title = "Biomasa vs. Zinc", x = "Zinc (ppm)", y = "Biomasa (g)")

(disp_pH | disp_sal) / (disp_zn | plot_spacer()) +
  plot_annotation(title = "Figura 2. Relación bivariada entre la biomasa y las covariables pH, Salinidad y Zinc",
                   theme = theme(plot.title = element_text(face = "bold", size = 12)))

Interpretación — Análisis bivariado Contrario a lo que podría anticiparse por el nombre de la base de datos, la covariable con mayor relación con la biomasa es el pH del suelo, con una correlación positiva muy fuerte y altamente significativa (r ≈ 0.93): a medida que el suelo se hace menos ácido (pH más cercano a la neutralidad), la producción de biomasa aumenta de forma marcada. Esto es biológicamente coherente, ya que la acidez extrema del suelo (en este conjunto de datos el pH oscila entre 3.2 y 7.45, es decir, de fuertemente ácido a neutro) reduce la disponibilidad de nutrientes esenciales, favorece la toxicidad por aluminio y manganeso, y limita la actividad microbiana del suelo, todo lo cual restringe severamente el crecimiento vegetal. El Zinc muestra también una relación fuerte pero negativa con la biomasa (r ≈ -0.78), un resultado que a primera vista podría parecer contraintuitivo (el zinc es un micronutriente esencial), pero que es consistente con un escenario de toxicidad por exceso de zinc en los suelos más ácidos de esta muestra: la acidez del suelo incrementa la solubilidad y biodisponibilidad de metales como el zinc, de modo que las muestras más ácidas (menor pH, menor biomasa) tienden a presentar simultáneamente mayor Zinc disponible, generando esta correlación negativa que refleja más un patrón de covarianza edáfica que un efecto tóxico directo aislado del zinc. La Salinidad, pese a dar nombre a la base de datos, no muestra relación apreciable con la biomasa en este conjunto de muestras (r ≈ -0.07, no significativa), y el Potasio tampoco exhibe una relación relevante (r ≈ -0.07, no significativa). En consecuencia, el pH se selecciona como la variable con mayor relación con la biomasa para el análisis de varianza del literal (c).

1.3 c) ANOVA de una vía: efecto del nivel de pH sobre la Biomasa

1.3.1 Categorización en terciles

terciles <- quantile(datos_salinidad$pH, probs = c(1/3, 2/3))

datos_salinidad <- datos_salinidad %>%
  mutate(pH_cat = case_when(
    pH <= terciles[1] ~ "Bajo",
    pH <= terciles[2] ~ "Medio",
    TRUE ~ "Alto"
  ) %>% factor(levels = c("Bajo", "Medio", "Alto")))

tibble(Tercil = c("1/3 (corte Bajo–Medio)", "2/3 (corte Medio–Alto)"),
       `Valor de pH` = round(terciles, 3)) %>%
  tabla_bonita("Tabla 3. Puntos de corte de terciles usados para categorizar el pH")
Tabla 3. Puntos de corte de terciles usados para categorizar el pH
Tercil Valor de pH
1/3 (corte Bajo–Medio) 3.883
2/3 (corte Medio–Alto) 4.900

Justificación del criterio de categorización: se utilizan los terciles empíricos (percentiles 33.3 % y 66.7 %) de la distribución de pH para definir los niveles Bajo, Medio y Alto. Este criterio es preferible a puntos de corte fijos arbitrarios (por ejemplo, umbrales agronómicos estándar de acidez) porque garantiza grupos de tamaño balanceado (n = 15 por grupo), lo cual maximiza la potencia estadística del ANOVA y evita la formación de grupos con muy pocas observaciones, a la vez que respeta la distribución real de los datos muestreados en este estudio particular.

resumen_grupos_p1 <- datos_salinidad %>%
  group_by(pH_cat) %>%
  summarise(n = n(),
            `Media Biomasa (g)` = round(mean(Biomasa), 2),
            `DE` = round(sd(Biomasa), 2),
            `EE` = round(sd(Biomasa)/sqrt(n), 2),
            .groups = "drop")
tabla_bonita(resumen_grupos_p1, "Tabla 4. Estadísticos de Biomasa por nivel categorizado de pH")
Tabla 4. Estadísticos de Biomasa por nivel categorizado de pH
pH_cat n Media Biomasa (g) DE EE
Bajo 15 593.02 164.83 42.56
Medio 15 1048.11 244.91 63.24
Alto 15 1605.39 547.60 141.39 
ggplot(datos_salinidad, aes(pH_cat, Biomasa, fill = pH_cat)) +
  geom_boxplot(alpha = 0.85, width = 0.55, outlier.color = "#B71C1C") +
  geom_jitter(width = 0.08, alpha = 0.45, size = 1.8, color = "#232323") +
  scale_fill_manual(values = c("Bajo" = "#7A4B3A", "Medio" = "#D9A441", "Alto" = "#A9C5A0")) +
  labs(title = "Figura 3. Biomasa según nivel categorizado de pH",
       subtitle = "Categorización por terciles (Bajo / Medio / Alto)",
       x = "Nivel de pH", y = "Biomasa (g)") +
  theme(legend.position = "none")

1.3.2 Verificación de supuestos del modelo

modelo_p1 <- aov(Biomasa ~ pH_cat, data = datos_salinidad)

shapiro_p1 <- shapiro.test(residuals(modelo_p1))
levene_p1  <- car::leveneTest(Biomasa ~ pH_cat, data = datos_salinidad)

tabla_bonita(tabla_supuestos(shapiro_p1, levene_p1),
             "Tabla 5. Verificación de supuestos del ANOVA (Punto 1)")
Tabla 5. Verificación de supuestos del ANOVA (Punto 1)
Prueba Estadístico Valor p Decisión
Shapiro-Wilk (normalidad de residuales) 0.973 0.375 No se rechaza H₀ → residuales compatibles con normalidad
Levene (homogeneidad de varianzas) 8.675 0.001 Se rechaza H₀ → varianzas heterogéneas 
qq_data <- tibble(res = residuals(modelo_p1))
p_qq <- ggplot(qq_data, aes(sample = res)) +
  stat_qq(color = paleta_principal[1], size = 1.8) +
  stat_qq_line(color = "#D9A441", linewidth = 0.9) +
  labs(title = "Q-Q plot de residuales", x = "Cuantiles teóricos", y = "Cuantiles muestrales")

p_res_fit <- tibble(fit = fitted(modelo_p1), res = residuals(modelo_p1)) %>%
  ggplot(aes(fit, res)) +
  geom_point(color = paleta_principal[2], size = 2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "#B71C1C") +
  labs(title = "Residuales vs. ajustados", x = "Valores ajustados", y = "Residuales")

p_qq | p_res_fit

Verificación de supuestos La prueba de Shapiro-Wilk sobre los residuales arroja un valor p = 0.3749, superior a 0.05, por lo que no se rechaza el supuesto de normalidad; el gráfico Q-Q confirma un ajuste razonable a la recta teórica. Sin embargo, la prueba de Levene arroja un valor p = 7^{-4}, claramente inferior a 0.05, por lo que se rechaza el supuesto de homogeneidad de varianzas: el grupo de pH “Alto” presenta una dispersión de biomasa considerablemente mayor que los grupos “Bajo” y “Medio” (véase la Tabla 4 y la Figura 3), lo cual es consistente con una mayor heterogeneidad biológica en las condiciones de suelo menos limitantes. Dado que no se cumple el supuesto de homogeneidad de varianzas, no es apropiado interpretar el ANOVA de Fisher clásico ni aplicar la prueba de LSD estándar, ya que ambos asumen varianzas iguales entre grupos. En su lugar, se emplea la prueba de Welch (ANOVA robusto a heterogeneidad de varianzas, que ajusta los grados de libertad mediante la corrección de Welch-Satterthwaite) junto con la prueba post-hoc de Games-Howell, que es el análogo robusto de la prueba de LSD para varianzas desiguales y tamaños de grupo posiblemente distintos.

1.3.3 Resultado del ANOVA (Welch, robusto a heterogeneidad de varianzas)

welch_p1 <- oneway.test(Biomasa ~ pH_cat, data = datos_salinidad, var.equal = FALSE)
tabla_welch_p1 <- tibble(
  `Estadístico F` = round(unname(welch_p1$statistic), 3),
  `gl numerador` = round(unname(welch_p1$parameter[1]), 2),
  `gl denominador` = round(unname(welch_p1$parameter[2]), 2),
  `Valor p` = signif(welch_p1$p.value, 4)
)
tabla_bonita(tabla_welch_p1, "Tabla 6. ANOVA de Welch — Biomasa ~ Nivel de pH (robusto a heterogeneidad de varianzas)")
Tabla 6. ANOVA de Welch — Biomasa ~ Nivel de pH (robusto a heterogeneidad de varianzas)
Estadístico F gl numerador gl denominador Valor p
34.727 2 24.74 0 
eta2_p1 <- effectsize::eta_squared(aov(Biomasa ~ pH_cat, data = datos_salinidad), partial = FALSE)
tabla_bonita(as.data.frame(eta2_p1), "Tabla 7. Tamaño del efecto (Eta cuadrado, calculado sobre el modelo ANOVA base)")
Tabla 7. Tamaño del efecto (Eta cuadrado, calculado sobre el modelo ANOVA base)
Parameter Eta2 CI CI_low CI_high
pH_cat 0.587 0.95 0.416 1

1.3.4 Comparaciones post-hoc (prueba de Games-Howell)

gh_p1 <- rstatix::games_howell_test(datos_salinidad, Biomasa ~ pH_cat)
tabla_bonita(gh_p1 %>% select(group1, group2, estimate, conf.low, conf.high, p.adj) %>%
               rename(`Grupo 1` = group1, `Grupo 2` = group2, `Diferencia de medias` = estimate,
                      `IC 95% inf.` = conf.low, `IC 95% sup.` = conf.high, `Valor p ajustado` = p.adj),
             "Tabla 8. Comparaciones post-hoc de Games-Howell entre niveles de pH (α = 0.05)")
Tabla 8. Comparaciones post-hoc de Games-Howell entre niveles de pH (α = 0.05)
Grupo 1 Grupo 2 Diferencia de medias IC 95% inf. IC 95% sup. Valor p ajustado
Bajo Medio 455.086 264.998 645.174 0.000
Bajo Alto 1012.362 632.537 1392.187 0.000
Medio Alto 557.276 164.440 950.112 0.005 
letras_p1 <- tibble(pH_cat = factor(c("Bajo","Medio","Alto"), levels = c("Bajo","Medio","Alto")),
                     letra = c("a", "b", "c"))

ggplot(datos_salinidad, aes(pH_cat, Biomasa, fill = pH_cat)) +
  geom_boxplot(alpha = 0.85, width = 0.5, outlier.shape = NA) +
  geom_jitter(width = 0.08, alpha = 0.35, size = 1.6) +
  geom_text(data = letras_p1, aes(y = max(datos_salinidad$Biomasa) * 1.08, label = letra),
            size = 5.5, fontface = "bold", color = "#1B4332") +
  scale_fill_manual(values = c("Bajo" = "#7A4B3A", "Medio" = "#D9A441", "Alto" = "#A9C5A0")) +
  labs(title = "Figura 4. Comparaciones post-hoc de Games-Howell entre niveles de pH",
       subtitle = "Letras distintas = diferencia estadísticamente significativa (α = 0.05, todos los pares difieren)",
       x = "Nivel de pH", y = "Biomasa (g)") +
  theme(legend.position = "none")

Conclusión — Punto 1. El ANOVA de Welch indica que el nivel categorizado de pH genera diferencias altamente significativas en la producción de biomasa de la planta forrajera (F = 34.73, p < 0.001), con un tamaño de efecto correspondiente a 58.7 % de la varianza de la biomasa asociada al nivel de pH. Las comparaciones post-hoc de Games-Howell muestran que los tres niveles difieren significativamente entre sí (todos los valores p ajustados < 0.01): las parcelas con pH Alto (suelo menos ácido, más cercano a la neutralidad) producen la mayor biomasa y con la mayor variabilidad, las de pH Medio una biomasa intermedia, y las de pH Bajo (suelo más ácido) la menor biomasa y la menor dispersión. Desde el punto de vista agronómico, esto indica que la acidez del suelo —y no su salinidad ni su contenido de potasio— es el principal factor limitante de la producción de esta especie forrajera dentro del rango de condiciones muestreadas, probablemente mediado por la reducción de disponibilidad de nutrientes y el incremento de la biodisponibilidad de metales potencialmente fitotóxicos (como el Zinc) en los suelos más ácidos. Es importante resaltar el valor pedagógico de este resultado: la variable que da nombre a la base de datos (Salinidad) no fue la de mayor relación con la biomasa, lo cual subraya la importancia de realizar siempre un análisis exploratorio bivariado objetivo antes de asumir cuál covariable es la más relevante.

2 Punto 2 — Consumo de oxígeno en moluscos bajo distintas concentraciones de agua de mar

Dos tipos de moluscos (A y B) fueron sometidos a tres concentraciones de agua de mar (50 %, 75 % y 100 %), registrándose el consumo de oxígeno (proporción de O₂ por unidad de peso seco), con 8 réplicas por combinación experimental (n = 48 en total).

load("moluscos.RData")
datos_moluscos <- BD_moluscos %>%
  mutate(Concentracion = factor(paste0(c_agua, "%"), levels = c("50%","75%","100%")),
         Tipo = factor(molusco, levels = c("A","B")),
         O2 = cons_o)

glimpse(datos_moluscos)
## Rows: 48
## Columns: 6
## $ c_agua        <dbl> 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 1…
## $ molusco       <chr> "A", "A", "A", "A", "A", "A", "A", "A", "B", "B", "B", "…
## $ cons_o        <dbl> 7.16, 8.26, 6.78, 14.00, 13.60, 11.10, 8.93, 9.66, 6.14,…
## $ Concentracion <fct> 100%, 100%, 100%, 100%, 100%, 100%, 100%, 100%, 100%, 10…
## $ Tipo          <fct> A, A, A, A, A, A, A, A, B, B, B, B, B, B, B, B, A, A, A,…
## $ O2            <dbl> 7.16, 8.26, 6.78, 14.00, 13.60, 11.10, 8.93, 9.66, 6.14,…

2.1 a) Análisis exploratorio univariado

resumen_p2_O2 <- resumen_descriptivo(datos_moluscos$O2) %>% mutate(Variable = "Consumo de O2 (global)", .before = 1)
tabla_bonita(resumen_p2_O2, "Tabla 9. Estadísticos descriptivos del consumo de oxígeno (muestra completa, n = 48)")
Tabla 9. Estadísticos descriptivos del consumo de oxígeno (muestra completa, n = 48)
Variable n Media DE CV (%) Mediana Mínimo Máximo Asimetría Curtosis
Consumo de O2 (global) 48 9.3 3.68 39.6 9.7 1.8 18.8 0.41 -0.19 
tabla_frecuencias_p2 <- datos_moluscos %>% count(Tipo, Concentracion, name = "n") %>%
  pivot_wider(names_from = Concentracion, values_from = n)
tabla_bonita(tabla_frecuencias_p2, "Tabla 10. Número de observaciones por combinación Tipo × Concentración (diseño balanceado, 8 réplicas por celda)")
Tabla 10. Número de observaciones por combinación Tipo × Concentración (diseño balanceado, 8 réplicas por celda)
Tipo 50% 75% 100%
A 8 8 8
B 8 8 8 
p2_hist <- ggplot(datos_moluscos, aes(O2)) +
  geom_histogram(aes(y = after_stat(density)), bins = 14, fill = paleta_principal[1], color = "white", alpha = 0.85) +
  geom_density(color = "#232323", linewidth = 0.7) +
  labs(title = "Distribución global del consumo de O2", x = "Consumo de O2 (proporción/peso seco)", y = "Densidad")

p2_box_tipo <- ggplot(datos_moluscos, aes(Tipo, O2, fill = Tipo)) +
  geom_boxplot(alpha = 0.85, width = 0.5) +
  scale_fill_manual(values = paleta_molusco) +
  labs(title = "Consumo de O2 por tipo de molusco", x = "Tipo", y = "Consumo de O2") +
  theme(legend.position = "none")

p2_hist | p2_box_tipo

Interpretación — Análisis univariado La distribución global del consumo de oxígeno muestra una dispersión considerable (rango aproximado 1.8–18.8), reflejo de que la muestra combina seis subgrupos experimentales (2 tipos × 3 concentraciones) con medias distintas; esto anticipa que buena parte de la variabilidad total se explica por los factores experimentales y no por dispersión aleatoria pura. El boxplot por tipo de molusco muestra medianas y dispersiones muy similares entre A y B, sin una diferencia visualmente marcada entre especies, lo que sugiere de entrada que el tipo de molusco podría no ser, por sí solo, un factor determinante del consumo de oxígeno — hipótesis que se contrasta formalmente con el ANOVA factorial en el literal (c).

2.2 b) Análisis exploratorio bivariado

resumen_bivariado_p2 <- datos_moluscos %>%
  group_by(Tipo, Concentracion) %>%
  summarise(n = n(), Media = round(mean(O2),2), DE = round(sd(O2),2), EE = round(sd(O2)/sqrt(n),3), .groups = "drop")
tabla_bonita(resumen_bivariado_p2, "Tabla 11. Consumo de O2 (media ± DE) por combinación Tipo × Concentración")
Tabla 11. Consumo de O2 (media ± DE) por combinación Tipo × Concentración
Tipo Concentracion n Media DE EE
A 50% 8 12.18 3.09 1.093
A 75% 8 7.89 2.74 0.969
A 100% 8 9.94 2.75 0.972
B 50% 8 12.33 3.52 1.244
B 75% 8 6.10 2.74 0.968
B 100% 8 7.41 2.84 1.006 
p2_box_conc <- ggplot(datos_moluscos, aes(Concentracion, O2, fill = Concentracion)) +
  geom_boxplot(alpha = 0.85, width = 0.55) +
  scale_fill_manual(values = c("50%"="#7A4B3A","75%"="#D9A441","100%"="#A9C5A0")) +
  labs(title = "Consumo de O2 según concentración de agua de mar (ambos tipos combinados)",
       x = "Concentración de agua de mar", y = "Consumo de O2") +
  theme(legend.position = "none")

p2_box_facet <- ggplot(datos_moluscos, aes(Concentracion, O2, fill = Tipo)) +
  geom_boxplot(alpha = 0.85, width = 0.6, position = position_dodge(0.7)) +
  scale_fill_manual(values = paleta_molusco) +
  labs(title = "Consumo de O2 según concentración, separado por tipo de molusco",
       x = "Concentración de agua de mar", y = "Consumo de O2")

interaccion_p2 <- datos_moluscos %>%
  group_by(Tipo, Concentracion) %>%
  summarise(Media = mean(O2), EE = sd(O2)/sqrt(n()), .groups = "drop")

p2_interaccion <- ggplot(interaccion_p2, aes(Concentracion, Media, color = Tipo, group = Tipo)) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 3.2) +
  geom_errorbar(aes(ymin = Media - EE, ymax = Media + EE), width = 0.08, linewidth = 0.7) +
  scale_color_manual(values = paleta_molusco) +
  labs(title = "Gráfico de interacción Tipo × Concentración",
       subtitle = "Líneas aproximadamente paralelas sugieren ausencia de interacción relevante",
       x = "Concentración de agua de mar", y = "Consumo medio de O2 (± EE)")

p2_box_conc / p2_box_facet / p2_interaccion +
  plot_annotation(title = "Figura 5. Análisis bivariado del consumo de oxígeno según concentración y tipo de molusco",
                   theme = theme(plot.title = element_text(face = "bold", size = 12)))

Interpretación — Análisis bivariado El consumo de oxígeno no sigue una relación monótona simple con la concentración de agua de mar: es más alto en la concentración 50 %, cae de forma marcada en 75 %, y se recupera parcialmente en 100 %, sin llegar al nivel observado en 50 %. Este patrón en forma de “valle” es fisiológicamente interpretable: en 50 % (medio fuertemente hiposalino) los moluscos probablemente activan mecanismos de compensación osmótica que elevan su tasa metabólica y, por tanto, su consumo de oxígeno; en 100 % (agua de mar sin diluir, su condición fisiológica de referencia) el consumo se mantiene en un nivel moderado correspondiente al metabolismo basal; mientras que en 75 % —una condición de estrés intermedio, posiblemente insuficiente para desencadenar una respuesta compensatoria plena pero suficiente para alterar la eficiencia metabólica— se observa el consumo más bajo de las tres condiciones. Este patrón es consistente entre los dos tipos de molusco (A y B siguen el mismo ordenamiento 50 % > 100 % > 75 %, aproximadamente), y el gráfico de interacción muestra líneas razonablemente paralelas entre A y B a través de las tres concentraciones, sin un cruce marcado ni una divergencia sistemática de pendientes. Esto sugiere, a nivel exploratorio, que la concentración de agua de mar afecta de forma similar a ambos tipos de molusco, sin evidencia visual fuerte de una interacción relevante — hipótesis que se contrasta formalmente en el literal (c).

2.3 c) ANOVA de dos vías: Tipo × Concentración

2.3.1 Verificación de supuestos del modelo

modelo_p2 <- aov(O2 ~ Tipo * Concentracion, data = datos_moluscos)

shapiro_p2 <- shapiro.test(residuals(modelo_p2))
levene_p2  <- car::leveneTest(O2 ~ Tipo * Concentracion, data = datos_moluscos)

tabla_bonita(tabla_supuestos(shapiro_p2, levene_p2),
             "Tabla 12. Verificación de supuestos del ANOVA de dos vías (Punto 2)")
Tabla 12. Verificación de supuestos del ANOVA de dos vías (Punto 2)
Prueba Estadístico Valor p Decisión
Shapiro-Wilk (normalidad de residuales) 0.958 0.086 No se rechaza H₀ → residuales compatibles con normalidad
Levene (homogeneidad de varianzas) 0.172 0.972 No se rechaza H₀ → varianzas homogéneas 
p2_qq <- tibble(res = residuals(modelo_p2)) %>%
  ggplot(aes(sample = res)) +
  stat_qq(color = paleta_principal[1], size = 1.8) +
  stat_qq_line(color = "#D9A441", linewidth = 0.9) +
  labs(title = "Q-Q plot de residuales", x = "Cuantiles teóricos", y = "Cuantiles muestrales")

p2_res_fit <- tibble(fit = fitted(modelo_p2), res = residuals(modelo_p2)) %>%
  ggplot(aes(fit, res)) +
  geom_point(color = paleta_principal[2], size = 2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "#B71C1C") +
  labs(title = "Residuales vs. ajustados", x = "Valores ajustados", y = "Residuales")

p2_qq | p2_res_fit

Verificación de supuestos La prueba de Shapiro-Wilk sobre los residuales del modelo factorial arroja p = 0.0857, por lo que no se rechaza la normalidad de los residuales, y la prueba de Levene arroja p = 0.9715, sin evidencia de heterogeneidad de varianzas entre las seis combinaciones experimentales. Ambos supuestos se satisfacen, por lo que el ANOVA factorial clásico y la prueba de LSD son procedimientos apropiados para este conjunto de datos.

2.3.2 Resultado del ANOVA factorial

tabla_anova_p2 <- broom::tidy(modelo_p2) %>% mutate(across(where(is.numeric), ~round(.x, 5)))
tabla_bonita(tabla_anova_p2, "Tabla 13. Tabla ANOVA de dos vías — Consumo de O2 ~ Tipo × Concentración")
Tabla 13. Tabla ANOVA de dos vías — Consumo de O2 ~ Tipo × Concentración
term df sumsq meansq statistic p.value
Tipo 1 23.227 23.227 2.651 0.111
Concentracion 2 230.816 115.408 13.171 0.000
Tipo:Concentracion 2 15.356 7.678 0.876 0.424
Residuals 42 368.011 8.762 NA NA 
eta2_p2 <- effectsize::eta_squared(modelo_p2, partial = TRUE)
tabla_bonita(as.data.frame(eta2_p2), "Tabla 14. Tamaños de efecto (Eta cuadrado parcial) por factor")
Tabla 14. Tamaños de efecto (Eta cuadrado parcial) por factor
Parameter Eta2_partial CI CI_low CI_high
Tipo 0.059 0.95 0.000 1
Concentracion 0.385 0.95 0.185 1
Tipo:Concentracion 0.040 0.95 0.000 1

Interpretación de los efectos principales y de la interacción El efecto principal del Tipo de molusco no resulta estadísticamente significativo (p = 0.111): con los datos disponibles, no hay evidencia suficiente de que el consumo medio de oxígeno difiera entre los moluscos tipo A y tipo B cuando se promedia sobre las tres concentraciones. El efecto principal de la Concentración de agua de mar, en cambio, es altamente significativo (p < 0.001) y es, con amplio margen, el factor que más varianza explica del consumo de oxígeno, confirmando el patrón en forma de “valle” (50 % > 100 % > 75 %) observado en el análisis bivariado. La interacción Tipo × Concentración tampoco resulta significativa (p = 0.4238), lo que confirma formalmente lo sugerido por el gráfico de interacción: el efecto de la concentración de agua de mar sobre el consumo de oxígeno es estadísticamente el mismo para ambos tipos de molusco. En consecuencia, el modelo se simplifica de forma natural a un análisis del efecto principal de la Concentración, sin necesidad de examinar las combinaciones específicas Tipo × Concentración por separado, y las comparaciones post-hoc se concentran en los tres niveles de concentración.

2.3.3 Comparaciones post-hoc (prueba de LSD) sobre el factor Concentración

Dado que ni el efecto principal de Tipo ni la interacción Tipo × Concentración resultaron significativos, las comparaciones post-hoc se concentran exclusivamente en el factor Concentración, que es el único efecto significativo del modelo factorial.

modelo_p2_conc <- aov(O2 ~ Concentracion, data = datos_moluscos)
lsd_p2 <- agricolae::LSD.test(modelo_p2_conc, "Concentracion", p.adj = "none", console = FALSE)

tabla_lsd_p2 <- lsd_p2$groups %>% rownames_to_column("Concentración de agua de mar") %>%
  rename(`Media O2` = O2, `Grupo LSD` = groups) %>%
  arrange(desc(`Media O2`))
tabla_bonita(tabla_lsd_p2, "Tabla 15. Comparaciones post-hoc LSD entre los tres niveles de Concentración de agua de mar")
Tabla 15. Comparaciones post-hoc LSD entre los tres niveles de Concentración de agua de mar
Concentración de agua de mar Media O2 Grupo LSD
50% 12.251 a
100% 8.671 b
75% 6.992 b 
grupos_lsd_p2 <- lsd_p2$groups %>% rownames_to_column("Concentracion") %>%
  mutate(Concentracion = factor(Concentracion, levels = c("50%","75%","100%")))

datos_moluscos %>%
  group_by(Concentracion) %>%
  summarise(Media = mean(O2), EE = sd(O2)/sqrt(n()), .groups = "drop") %>%
  left_join(grupos_lsd_p2, by = "Concentracion") %>%
  ggplot(aes(Concentracion, Media, fill = Concentracion)) +
  geom_col(alpha = 0.88, width = 0.55) +
  geom_errorbar(aes(ymin = Media - EE, ymax = Media + EE), width = 0.15) +
  geom_text(aes(label = groups, y = Media + EE + 0.6), size = 5.5, fontface = "bold", color = "#1B4332") +
  scale_fill_manual(values = c("50%"="#7A4B3A","75%"="#D9A441","100%"="#A9C5A0")) +
  labs(title = "Figura 6. Comparaciones post-hoc LSD entre niveles de Concentración de agua de mar",
       subtitle = "Letras distintas = diferencia estadísticamente significativa (α = 0.05); ambos tipos de molusco combinados",
       x = "Concentración de agua de mar", y = "Consumo medio de O2 (± EE)") +
  theme(legend.position = "none")

Conclusión — Punto 2. El ANOVA factorial muestra que solo la concentración de agua de mar afecta significativamente el consumo de oxígeno de los moluscos; ni el tipo de molusco por sí solo, ni la interacción entre tipo y concentración, alcanzan significancia estadística con los datos disponibles. Las comparaciones post-hoc LSD sobre el factor Concentración indican que el consumo de oxígeno en 50 % es significativamente mayor que en 75 % y que en 100 %, mientras que 75 % y 100 % no difieren significativamente entre sí. Biológicamente, este patrón sugiere que la dilución fuerte del agua de mar (50 %) impone un estrés osmótico suficientemente severo como para desencadenar un incremento metabólico compensatorio medible en ambas especies por igual, mientras que la reducción moderada de la concentración (75 %) no genera una respuesta compensatoria clara y podría incluso reflejar una fase de menor eficiencia metabólica transitoria antes de que se active la compensación plena. La ausencia de diferencias entre tipos de molusco indica que, en el rango de concentraciones evaluado, A y B comparten una capacidad osmorregulatoria estadísticamente equivalente frente al estrés hiposalino, un resultado relevante para entender su tolerancia comparativa a la variabilidad salina de su hábitat.

3 Punto 3 — Efecto del uso del suelo sobre la biodiversidad de anfibios

Se evalúa el efecto del uso del suelo sobre la biodiversidad de anfibios en cuatro tipos de hábitat (bosque primario, bosque secundario, sistema silvopastoril y potrero), representando un gradiente de intervención antrópica, con 52 parcelas de muestreo (13 por hábitat).

load("Biodiversidad.RData")
datos_biodiversidad <- BD_biodiversidad %>%
  mutate(Habitat = factor(Habitat, levels = c("Bosque primario","Bosque secundario","Sistema silvopastoril","Potrero")))

glimpse(datos_biodiversidad)
## Rows: 52
## Columns: 5
## $ Parcela <chr> "P001", "P002", "P003", "P004", "P005", "P006", "P007", "P008"…
## $ Habitat <fct> Bosque primario, Bosque primario, Bosque primario, Bosque prim…
## $ Riqueza <dbl> 14, 14, 13, 14, 10, 15, 19, 16, 15, 21, 14, 13, 14, 10, 14, 15…
## $ Shannon <dbl> 2.13, 2.98, 2.07, 2.53, 2.82, 2.75, 2.45, 1.73, 2.04, 2.48, 2.…
## $ Altitud <dbl> 1478.0, 1476.8, 1505.2, 1395.2, 1493.5, 1337.9, 1482.4, 1600.8…

3.1 a) Análisis exploratorio univariado

resumen_p3 <- bind_rows(
  resumen_descriptivo(datos_biodiversidad$Riqueza) %>% mutate(Variable = "Riqueza de especies", .before = 1),
  resumen_descriptivo(datos_biodiversidad$Shannon)  %>% mutate(Variable = "Índice de Shannon-Wiener", .before = 1),
  resumen_descriptivo(datos_biodiversidad$Altitud)  %>% mutate(Variable = "Altitud (m.s.n.m.)", .before = 1)
)
tabla_bonita(resumen_p3, "Tabla 16. Estadísticos descriptivos globales (n = 52 parcelas)")
Tabla 16. Estadísticos descriptivos globales (n = 52 parcelas)
Variable n Media DE CV (%) Mediana Mínimo Máximo Asimetría Curtosis
Riqueza de especies 52 9.73 4.37 44.9 9.50 2.00 21.00 0.29 -0.63
Índice de Shannon-Wiener 52 1.76 0.61 34.4 1.85 0.67 2.98 0.04 -0.75
Altitud (m.s.n.m.) 52 1243.89 172.71 13.9 1254.80 950.90 1600.80 0.09 -1.19 
p3_hist_riq <- ggplot(datos_biodiversidad, aes(Riqueza)) +
  geom_histogram(aes(y = after_stat(density)), bins = 12, fill = paleta_principal[1], color = "white", alpha = 0.85) +
  geom_density(color = "#232323", linewidth = 0.7) +
  labs(title = "Riqueza de especies", x = "N° de especies", y = "Densidad")

p3_hist_shan <- ggplot(datos_biodiversidad, aes(Shannon)) +
  geom_histogram(aes(y = after_stat(density)), bins = 12, fill = paleta_principal[4], color = "white", alpha = 0.85) +
  geom_density(color = "#232323", linewidth = 0.7) +
  labs(title = "Índice de Shannon-Wiener", x = "H'", y = "Densidad")

p3_hist_alt <- ggplot(datos_biodiversidad, aes(Altitud)) +
  geom_histogram(aes(y = after_stat(density)), bins = 12, fill = paleta_principal[2], color = "white", alpha = 0.85) +
  geom_density(color = "#232323", linewidth = 0.7) +
  labs(title = "Altitud", x = "m.s.n.m.", y = "Densidad")

p3_box_global <- datos_biodiversidad %>%
  select(Riqueza, Shannon) %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "Valor") %>%
  ggplot(aes(Variable, Valor, fill = Variable)) +
  geom_boxplot(alpha = 0.85, width = 0.5) +
  scale_fill_manual(values = c("Riqueza" = paleta_principal[1], "Shannon" = paleta_principal[4])) +
  labs(title = "Boxplots globales (Riqueza y Shannon)", x = NULL, y = "Valor") +
  theme(legend.position = "none")

(p3_hist_riq | p3_hist_shan) / (p3_hist_alt | p3_box_global) +
  plot_annotation(title = "Figura 7. Distribuciones univariadas de riqueza, diversidad de Shannon y altitud",
                   theme = theme(plot.title = element_text(face = "bold", size = 12)))

Interpretación — Análisis univariado La riqueza de especies muestra una distribución con una dispersión amplia (desde parcelas muy pobres hasta parcelas con alta riqueza), reflejando la heterogeneidad esperada entre los cuatro tipos de hábitat con distinto grado de intervención. El índice de Shannon-Wiener presenta una distribución razonablemente simétrica con dos concentraciones de valores (más evidente cuando se desagrega por hábitat, ver literal b), consistente con la coexistencia de ambientes muy diversos (bosques) y ambientes empobrecidos (potreros) dentro de la misma muestra global. La altitud se distribuye de forma aproximadamente normal en el rango de la reserva forestal muestreada, sin valores atípicos evidentes, lo que sugiere un gradiente altitudinal continuo y bien representado entre las 52 parcelas.

3.2 b) Análisis exploratorio bivariado

resumen_bivariado_p3 <- datos_biodiversidad %>%
  group_by(Habitat) %>%
  summarise(n = n(),
            `Riqueza media` = round(mean(Riqueza), 1), `DE Riqueza` = round(sd(Riqueza), 1),
            `Shannon medio` = round(mean(Shannon), 2), `DE Shannon` = round(sd(Shannon), 2),
            `Altitud media` = round(mean(Altitud), 0),
            .groups = "drop")
tabla_bonita(resumen_bivariado_p3, "Tabla 17. Riqueza, diversidad de Shannon y altitud por tipo de hábitat")
Tabla 17. Riqueza, diversidad de Shannon y altitud por tipo de hábitat
Habitat n Riqueza media DE Riqueza Shannon medio DE Shannon Altitud media
Bosque primario 13 14.8 2.7 2.47 0.38 1453
Bosque secundario 13 11.1 3.0 1.97 0.22 1326
Sistema silvopastoril 13 8.1 2.5 1.60 0.28 1151
Potrero 13 5.0 1.3 1.00 0.24 1045 
p3_box_riq_hab <- ggplot(datos_biodiversidad, aes(Habitat, Riqueza, fill = Habitat)) +
  geom_boxplot(alpha = 0.88, width = 0.6) +
  geom_jitter(width = 0.08, alpha = 0.4, size = 1.6) +
  scale_fill_manual(values = paleta_habitat) +
  labs(title = "Riqueza de especies por tipo de hábitat", x = NULL, y = "N° de especies") +
  theme(legend.position = "none", axis.text.x = element_text(angle = 18, hjust = 1))

p3_box_shan_hab <- ggplot(datos_biodiversidad, aes(Habitat, Shannon, fill = Habitat)) +
  geom_boxplot(alpha = 0.88, width = 0.6) +
  geom_jitter(width = 0.08, alpha = 0.4, size = 1.6) +
  scale_fill_manual(values = paleta_habitat) +
  labs(title = "Índice de Shannon-Wiener por tipo de hábitat", x = NULL, y = "H'") +
  theme(legend.position = "none", axis.text.x = element_text(angle = 18, hjust = 1))

p3_disp_alt_shan <- ggplot(datos_biodiversidad, aes(Altitud, Shannon, color = Habitat)) +
  geom_point(size = 2.6, alpha = 0.8) +
  geom_smooth(aes(group = 1), method = "lm", se = TRUE, color = "#232323", linewidth = 0.8, linetype = "dashed") +
  scale_color_manual(values = paleta_habitat) +
  labs(title = "Relación entre altitud y diversidad de Shannon",
       subtitle = "Línea de tendencia general (independiente del hábitat) en negro punteado",
       x = "Altitud (m.s.n.m.)", y = "Índice de Shannon-Wiener")

p3_box_riq_hab / p3_box_shan_hab / p3_disp_alt_shan +
  plot_annotation(title = "Figura 8. Análisis bivariado: riqueza y diversidad por hábitat, y relación altitud-diversidad",
                   theme = theme(plot.title = element_text(face = "bold", size = 12)))

cor_alt_shannon <- cor.test(datos_biodiversidad$Altitud, datos_biodiversidad$Shannon)
tibble(`r de Pearson` = round(unname(cor_alt_shannon$estimate), 3),
       `Valor p` = round(cor_alt_shannon$p.value, 5),
       IC_95 = paste0("[", round(cor_alt_shannon$conf.int[1],3), ", ", round(cor_alt_shannon$conf.int[2],3), "]")) %>%
  tabla_bonita("Tabla 18. Correlación entre altitud y diversidad de Shannon")
Tabla 18. Correlación entre altitud y diversidad de Shannon
r de Pearson Valor p IC_95
0.752 0 [0.603, 0.851]

Interpretación — Análisis bivariado Tanto la riqueza de especies como el índice de Shannon-Wiener muestran un patrón decreciente claro y consistente a lo largo del gradiente de intervención antrópica: bosque primario > bosque secundario > sistema silvopastoril > potrero. Este patrón es coherente con la teoría ecológica de que la conversión y simplificación estructural del hábitat (pérdida de estratos vegetales, de microhábitats húmedos y de conectividad) reduce tanto el número de especies de anfibios presentes como la equidad en su distribución. Respecto a la relación entre altitud y diversidad, se observa una correlación positiva fuerte y estadísticamente significativa (r ≈ 0.75): las parcelas ubicadas a mayor altitud dentro de la reserva tienden a presentar mayor diversidad de Shannon. Es importante señalar que esta asociación está en buena parte confundida con el tipo de hábitat, dado que en este diseño los bosques (más diversos) se ubican preferentemente en las zonas más altas y menos accesibles de la reserva, mientras que los potreros (menos diversos) ocupan las zonas bajas más intervenidas; por tanto, la altitud actúa aquí como un proxy parcial del gradiente de intervención antrópica más que como un factor ecológico independiente.

3.3 c) ANOVA de una vía: efecto del hábitat sobre el índice de Shannon-Wiener

3.3.1 Verificación de supuestos del modelo

modelo_p3 <- aov(Shannon ~ Habitat, data = datos_biodiversidad)

shapiro_p3 <- shapiro.test(residuals(modelo_p3))
levene_p3  <- car::leveneTest(Shannon ~ Habitat, data = datos_biodiversidad)

tabla_bonita(tabla_supuestos(shapiro_p3, levene_p3),
             "Tabla 19. Verificación de supuestos del ANOVA (Punto 3)")
Tabla 19. Verificación de supuestos del ANOVA (Punto 3)
Prueba Estadístico Valor p Decisión
Shapiro-Wilk (normalidad de residuales) 0.979 0.462 No se rechaza H₀ → residuales compatibles con normalidad
Levene (homogeneidad de varianzas) 1.277 0.293 No se rechaza H₀ → varianzas homogéneas 
p3_qq <- tibble(res = residuals(modelo_p3)) %>%
  ggplot(aes(sample = res)) +
  stat_qq(color = paleta_principal[1], size = 1.8) +
  stat_qq_line(color = "#D9A441", linewidth = 0.9) +
  labs(title = "Q-Q plot de residuales", x = "Cuantiles teóricos", y = "Cuantiles muestrales")

p3_res_fit <- tibble(fit = fitted(modelo_p3), res = residuals(modelo_p3)) %>%
  ggplot(aes(fit, res)) +
  geom_point(color = paleta_principal[2], size = 2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "#B71C1C") +
  labs(title = "Residuales vs. ajustados", x = "Valores ajustados", y = "Residuales")

p3_qq | p3_res_fit

Verificación de supuestos La prueba de Shapiro-Wilk sobre los residuales arroja p = 0.4625, sin evidencia en contra de la normalidad, y la prueba de Levene arroja p = 0.293, sin evidencia de heterogeneidad de varianzas entre hábitats. Ambos supuestos se cumplen satisfactoriamente, validando el uso del ANOVA de una vía clásico y de la prueba de LSD para las comparaciones post-hoc entre los cuatro tipos de hábitat.

3.3.2 Resultado del ANOVA

tabla_anova_p3 <- broom::tidy(modelo_p3) %>% mutate(across(where(is.numeric), ~round(.x, 6)))
tabla_bonita(tabla_anova_p3, "Tabla 20. Tabla ANOVA — Índice de Shannon-Wiener ~ Hábitat")
Tabla 20. Tabla ANOVA — Índice de Shannon-Wiener ~ Hábitat
term df sumsq meansq statistic p.value
Habitat 3 14.862 4.954 60.611 0
Residuals 48 3.923 0.082 NA NA 
eta2_p3 <- effectsize::eta_squared(modelo_p3, partial = FALSE)
tabla_bonita(as.data.frame(eta2_p3), "Tabla 21. Tamaño del efecto (Eta cuadrado)")
Tabla 21. Tamaño del efecto (Eta cuadrado)
Parameter Eta2 CI CI_low CI_high
Habitat 0.791 0.95 0.699 1

3.3.3 Comparaciones post-hoc (prueba de LSD)

lsd_p3 <- agricolae::LSD.test(modelo_p3, "Habitat", p.adj = "none", console = FALSE)
tabla_lsd_p3 <- lsd_p3$groups %>% rownames_to_column("Hábitat") %>%
  rename(`Shannon medio` = Shannon, `Grupo LSD` = groups)
tabla_bonita(tabla_lsd_p3, "Tabla 22. Comparaciones post-hoc LSD entre los cuatro tipos de hábitat")
Tabla 22. Comparaciones post-hoc LSD entre los cuatro tipos de hábitat
Hábitat Shannon medio Grupo LSD
Bosque primario 2.468 a
Bosque secundario 1.975 b
Sistema silvopastoril 1.603 c
Potrero 1.004 d 
grupos_lsd_p3 <- lsd_p3$groups %>% rownames_to_column("Habitat") %>%
  mutate(Habitat = factor(Habitat, levels = c("Bosque primario","Bosque secundario","Sistema silvopastoril","Potrero")))

ggplot(datos_biodiversidad, aes(Habitat, Shannon, fill = Habitat)) +
  geom_boxplot(alpha = 0.88, width = 0.55, outlier.shape = NA) +
  geom_jitter(width = 0.08, alpha = 0.35, size = 1.6) +
  geom_text(data = grupos_lsd_p3, aes(y = max(datos_biodiversidad$Shannon) * 1.12, label = groups),
            size = 5.5, fontface = "bold", color = "#1B4332") +
  scale_fill_manual(values = paleta_habitat) +
  labs(title = "Figura 9. Comparaciones post-hoc LSD del índice de Shannon-Wiener entre hábitats",
       subtitle = "Letras distintas = diferencia estadísticamente significativa (α = 0.05)",
       x = NULL, y = "Índice de Shannon-Wiener") +
  theme(legend.position = "none", axis.text.x = element_text(angle = 15, hjust = 1))

Conclusión — Punto 3. El ANOVA de una vía muestra que el tipo de hábitat afecta de manera altamente significativa el índice de diversidad de Shannon-Wiener de anfibios (F(3, 48) = 60.61, p < 0.001), explicando 79.1 % de la variabilidad total observada. Las comparaciones post-hoc LSD confirman que todos o casi todos los hábitats difieren significativamente entre sí, siguiendo el gradiente esperado de intervención antrópica: el bosque primario presenta la mayor diversidad, seguido del bosque secundario, luego el sistema silvopastoril, y finalmente el potrero con la menor diversidad. Este resultado aporta evidencia estadística robusta de que la transformación del uso del suelo —incluso en sistemas de manejo intermedio como el silvopastoril— tiene un costo medible sobre la biodiversidad de anfibios, un grupo particularmente sensible a la pérdida de humedad, cobertura vegetal y heterogeneidad estructural del hábitat. Estos hallazgos respaldan la relevancia de conservar remanentes de bosque primario y promover la restauración o el manejo silvopastoril como estrategias de mitigación frente a la conversión total a potrero.

Conclusiones generales del taller

Los tres análisis desarrollados en este informe ilustran, desde sistemas biológicos distintos, la utilidad del análisis de varianza como herramienta para contrastar hipótesis sobre el efecto de factores ambientales o experimentales categóricos sobre una variable respuesta continua, siempre que se verifiquen adecuadamente sus supuestos:

  1. En el Punto 1, el pH del suelo —y no la Salinidad, pese a dar nombre a la base de datos— emergió como el principal determinante de la producción de biomasa vegetal; al no cumplirse el supuesto de homogeneidad de varianzas, fue necesario recurrir al ANOVA de Welch y a la prueba post-hoc de Games-Howell, ambos robustos a la heterogeneidad de varianzas, confirmando diferencias significativas entre los tres niveles categorizados de pH.
  2. En el Punto 2, el diseño factorial mostró que únicamente la concentración de agua de mar afecta significativamente el consumo de oxígeno —siguiendo un patrón no monótono con mayor consumo bajo dilución fuerte (50 %)—, mientras que ni el tipo de molusco ni la interacción Tipo × Concentración resultaron significativos, indicando una respuesta osmorregulatoria estadísticamente equivalente entre las dos especies evaluadas.
  3. En el Punto 3, el gradiente de intervención antrópica sobre el uso del suelo mostró un efecto negativo, gradual y estadísticamente robusto sobre la diversidad de anfibios, reforzando la relevancia de la conservación de bosques primarios y del manejo silvopastoril como estrategias de conservación de la biodiversidad.

En los tres casos, la verificación previa de los supuestos de normalidad de residuales (Shapiro-Wilk) y homogeneidad de varianzas (Levene) resultó decisiva para elegir el procedimiento estadístico correcto: el ANOVA paramétrico clásico con LSD de Fisher cuando los supuestos se cumplieron (Puntos 2 y 3), y su alternativa robusta —Welch con Games-Howell— cuando la homogeneidad de varianzas no se sostuvo (Punto 1), garantizando en todos los casos la validez estadística de las conclusiones biológicas presentadas.

Informe elaborado en R Markdown · Universidad del Valle, Facultad de Ciencias Naturales y Exactas · Cali, Colombia