1 Introducción

Este informe recoge el desarrollo del Taller 1 del curso de Diseño de Experimentos. El taller está compuesto por tres problemas independientes, cada uno asociado a una base de datos propia (Salinidad.RData, moluscos.RData y Biodiversidad.RData), en los que se combina un análisis exploratorio de datos (AED) con un Análisis de Varianza (ANOVA) y, cuando la prueba global resulta significativa, se realizan análisis complementarios usando pruebas de comparación múltiple post-ANOVA según el requerimiento.

En términos generales, para cada uno de los tres puntos se sigue la misma ruta de trabajo:

  1. Exploración univariada: estadísticos descriptivos (tendencia central, dispersión, forma) y gráficos (histogramas, densidades, boxplots) para cada variable involucrada.
  2. Exploración bivariada: relación entre la variable respuesta y las covariables/factores de interés, mediante diagramas de dispersión, correlaciones o boxplots comparativos.
  3. Modelamiento vía ANOVA: planteamiento del modelo (una vía o dos vías, según el diseño), con su respectiva tabla de análisis de varianza.
  4. Verificación de supuestos: normalidad de los residuales (prueba de Shapiro-Wilk, complementada con un gráfico Q-Q) y homogeneidad de varianzas (prueba de Levene), todos los 3 conjuntos de datos se representan mediante el modelo de medias.
  5. Comparaciones múltiples: cuando el ANOVA resulta significativo, se aplica la prueba correspondiente para identificar específicamente cuáles grupos difieren entre sí.
  6. Conclusiones puntuales.

Todos los análisis se realizaron con R (funciones base aov, lm, shapiro.test, pairwise.t.test) y los paquetes ggplot2, car, dplyr, GGally y multcompView, fijando un nivel de significancia \(\alpha = 0.05\) en todas las pruebas de hipótesis.

2 Punto 1: Biomasa forrajera y características del suelo

2.1 Introducción

La producción de biomasa de una planta forrajera está relacionada con las condiciones ambientales en las que se desarrolla, entre ellas las características fisicoquímicas del suelo. Factores como el pH, la salinidad y la disponibilidad de nutrientes y micronutrientes pueden afectar procesos fundamentales para el crecimiento vegetal, como la absorción de agua y minerales, la actividad enzimática y el desarrollo de los tejidos (Neina, 2019). Por esta razón, para conocer si la producción de biomasa de la planta forrajera cambia dependiendo de las condiciones del suelo, e identificar cual es el factor que más influye se realizó el siguiente análisis estadístico.

2.2 Metodología

Se dispone de 45 muestras de suelo en las que se midió la producción de biomasa de una planta forrajera (variable respuesta, en gramos) y cuatro covariables del suelo en el que crecía: pH, Salinidad, Zinc y Potasio. El objetivo es (a) describir cada variable, (b) identificar cuál covariable se relaciona más fuertemente con la biomasa y (c) categorizar dicha covariable en tres niveles (bajo/medio/alto, mediante terciles) para evaluarla como factor en un ANOVA de una vía, con LSD de Fisher como prueba post-hoc.

load("Salinidad.RData")
d1 <- Salinidad
kable(head(d1, 6), caption = "Primeras filas de la base Salinidad")
Primeras filas de la base Salinidad
Biomasa pH Salinidad Zinc Potasio
765.280 5.00 33 16.4524 1441.67
954.017 4.70 35 13.9852 1299.19
827.686 4.20 32 15.3276 1154.27
755.072 4.40 30 17.3128 1045.15
896.176 5.55 33 22.3312 521.62
1422.836 5.50 33 12.2778 1273.02

2.3 Exploración de datos — análisis univariado

tabla_u1 <- rbind(
  resumen_num(d1$Biomasa,   "Biomasa (g)"),
  resumen_num(d1$pH,        "pH"),
  resumen_num(d1$Salinidad, "Salinidad"),
  resumen_num(d1$Zinc,      "Zinc"),
  resumen_num(d1$Potasio,   "Potasio")
)
kable(tabla_u1, caption = "Estadísticos descriptivos de la biomasa y sus covariables (n = 45)")
Estadísticos descriptivos de la biomasa y sus covariables (n = 45)
Variable n Media DE CV_pct Mínimo Q1 Mediana Q3 Máximo Asimetría
25% Biomasa (g) 45 1082.17 546.29 50.5 369.82 654.82 991.83 1346.88 2337.33 0.91
25%1 pH 45 4.61 1.25 27.2 3.20 3.45 4.45 5.35 7.45 0.87
25%2 Salinidad 45 30.27 3.72 12.3 24.00 27.00 30.00 33.00 38.00 0.31
25%3 Zinc 45 17.83 8.27 46.4 0.21 13.99 19.24 22.68 31.29 -0.66
25%4 Potasio 45 797.38 297.58 37.3 350.73 526.97 773.30 954.11 1441.67 0.48 
col5 <- rep(paleta, length.out = 5)
wrap_plots(
  poligono_box(d1$Biomasa,   "Biomasa (g)", col5[1]),
  poligono_box(d1$pH,        "pH",          col5[2]),
  poligono_box(d1$Salinidad, "Salinidad",   col5[3]),
  poligono_box(d1$Zinc,      "Zinc",        col5[4]),
  poligono_box(d1$Potasio,   "Potasio",     col5[5]),
  ncol = 3
) + plot_annotation(title = "Distribución univariada de la biomasa y las covariables del suelo",
                     theme = theme(plot.title = element_text(face="bold", color="#3B5D42", family="serif")))

La biomasa presenta una media de 1082,2g (DE=546,3) y un rango entre 369,8 y 2337,3 g, con un coeficiente de variación del 50,5%, lo que evidencia una dispersión considerable entre sitios y una asimetría positiva moderada (0,91). Si se comparan, se muestra que la mayoría de las parcelas producen valoress intermedios, y solo unas pocas alcanzan producciones altas que jalan la media hacia arriba. El pH del suelo osciló entre valores fuertemente ácidos (3.2) y cercanos a la neutralidad (7.45), con una distribución razonablemente simétrica, por lo que se asume que predominan los suelos ácidos, con algunos sitios con pH cercanos a la neutralidad. La Salinidad en cambio, se concentra en un rango estrecho (24-38) y con baja dispersión relativa, mientras que el Zinc y el Potasio muestran mayor variabilidad entre muestras, reflejo de la heterogeneidad ambiental de los sitios muestreados.

2.4 Exploración de datos — análisis bivariado

ggpairs(d1[, c("Biomasa","pH","Salinidad","Zinc")],
        lower = list(continuous = wrap("smooth", alpha = 0.5, color = paleta[1])),
        upper = list(continuous = wrap("cor", size = 4))) +
  labs(title = "Matriz de dispersión y correlaciones: Biomasa vs. covariables") +
  theme(strip.text = element_text(size = 9))

cor_bio <- cor(d1[, c("Biomasa","pH","Salinidad","Zinc")])[1, -1]
kable(data.frame(Covariable = names(cor_bio),
                  Correlación_con_Biomasa = round(cor_bio, 3)),
      caption = "Correlación de Pearson de cada covariable con la Biomasa")
Correlación de Pearson de cada covariable con la Biomasa
Covariable Correlación_con_Biomasa
pH pH 0.928
Salinidad Salinidad -0.067
Zinc Zinc -0.781

La covariable con mayor relación (en valor absoluto) con la biomasa es el pH (\(r =\) 0.928), seguida por el Zinc (\(r =\) -0.781); ambas relaciones son de magnitud considerable, pero de signo opuesto: a mayor pH (suelos menos ácidos) la biomasa tiende a aumentar, mientras que a mayor concentración de Zinc tiende a disminuir (posible efecto de fitotoxicidad). En contraste, resulta llamativo que la variable que da nombre a la base de datos, la Salinidad, es la que presenta menor relación lineal con la biomasa (\(r =\) -0.067), prácticamente nula, lo que también puede relacionarse con la baja dispersión de esta variable. Por lo anterior, la covariable seleccionada para construir el factor del ANOVA es el pH.

2.5 Categorización de la covariable y Análisis de Varianza

terciles <- quantile(d1$pH, probs = c(1/3, 2/3))
d1$pH_cat <- cut(d1$pH,
                  breaks = c(-Inf, terciles[1], terciles[2], Inf),
                  labels = c("Bajo","Medio","Alto"))
kable(as.data.frame(table(d1$pH_cat)) %>% setNames(c("Nivel de pH","n")),
      caption = paste0("Niveles de pH construidos por terciles (cortes en ",
                        round(terciles[1],2), " y ", round(terciles[2],2), ")"))
Niveles de pH construidos por terciles (cortes en 3.88 y 4.9)
Nivel de pH n
Bajo 15
Medio 15
Alto 15

Se categorizó el pH en tres niveles (bajo/medio/alto) usando los terciles muestrales (percentiles 33.3% y 66.7%) como puntos de corte, lo que garantiza grupos similares de 15 observaciones cada uno y evita la arbitrariedad de límites fijos no basados en los datos. En los boxplot de biomasa por nivel de pH se ve una separación clara de las cajas a lo largo de los tres niveles del gradiente de acidez, con la mediana de biomasa creciendo de forma constante de bajo a alto y un solapamiento mínimo entre las cajas.

ggplot(d1, aes(x = pH_cat, y = Biomasa, fill = pH_cat)) +
  geom_boxplot(alpha = 0.85, outlier.color = "#B3703A") +
  geom_jitter(width = 0.08, alpha = 0.5, size = 1.6, color = "#3A3833") +
  scale_fill_manual(values = paleta[c(3,5,1)]) +
  labs(title = "Biomasa según nivel categorizado de pH del suelo",
       x = "Nivel de pH", y = "Biomasa (g)") +
  theme(legend.position = "none")

2.5.1 Modelo ANOVA de una vía

\[Biomasa_{ij} = \mu + \tau_i + \varepsilon_{ij}, \qquad i = \text{Bajo, Medio, Alto}\]

mod1 <- aov(Biomasa ~ pH_cat, data = d1)
tabla_anova1 <- summary(mod1)[[1]]
kable(round(tabla_anova1, 4), caption = "Tabla ANOVA — Biomasa según nivel de pH")
Tabla ANOVA — Biomasa según nivel de pH
Df Sum Sq Mean Sq F value Pr(>F)
pH_cat 2 7712683 3856341.6 29.8928 0
Residuals 42 5418235 129005.6 NA NA 
eta2_1 <- tabla_anova1["pH_cat","Sum Sq"] / sum(tabla_anova1[,"Sum Sq"])

Resultado del ANOVA. El estadístico de prueba es \(F\)(2, 42) = 29.89 con un valor \(p\) = <1e-04, muy inferior a \(\alpha = 0.05\). Por lo tanto, se rechaza la hipótesis nula de igualdad de medias: el nivel de pH del suelo sí genera diferencias significativas en la producción de biomasa. El tamaño del efecto (\(\eta^2 =\) 0.587) indica que el nivel de pH explica cerca del 59% de la variabilidad total de la biomasa, un efecto grande.

2.6 Verificación de supuestos

par_grid <- gridExtra::arrangeGrob(
  ggplot(data.frame(res = residuals(mod1)), aes(sample = res)) +
    stat_qq(color = paleta[1]) + stat_qq_line(color = "#B3703A") +
    labs(title = "Q-Q plot de los residuales", x="Cuantiles teóricos", y="Cuantiles muestrales"),
  ggplot(data.frame(res = residuals(mod1), fit = fitted(mod1)), aes(x = fit, y = res)) +
    geom_point(color = paleta[1]) + geom_hline(yintercept = 0, linetype = "dashed", color="#B3703A") +
    labs(title = "Residuales vs. valores ajustados", x = "Valores ajustados", y = "Residuales"),
  ncol = 2
)
grid::grid.draw(par_grid)

sh1 <- shapiro.test(residuals(mod1))
lv1 <- leveneTest(Biomasa ~ pH_cat, data = d1)
kable(data.frame(Prueba = c("Shapiro-Wilk (normalidad residuales)", "Levene (homogeneidad de varianzas)"),
                  Estadístico = c(round(sh1$statistic,3), round(lv1$`F value`[1],3)),
                  Valor_p = c(format.pval(sh1$p.value, digits=3),
                              format.pval(lv1$`Pr(>F)`[1], digits=3))),
      caption = "Verificación de supuestos del modelo ANOVA — Punto 1")
Verificación de supuestos del modelo ANOVA — Punto 1
Prueba Estadístico Valor_p
W Shapiro-Wilk (normalidad residuales) 0.973 0.375
Levene (homogeneidad de varianzas) 8.675 0.000702
  • Normalidad: la prueba de Shapiro-Wilk sobre los residuales del modelo da \(p =\) 0.375, > 0.05, por lo que no se rechaza el supuesto de normalidad; el gráfico Q-Q corrobora un ajuste razonable a la recta teórica, sin colas marcadamente pesadas.
  • Homogeneidad de varianzas: la prueba de Levene arroja \(p =\) 0.000702, < 0.05, por lo que se rechaza el supuesto de homocedasticidad: la variabilidad de la biomasa aumenta con el nivel de pH (grupo “Alto” notablemente más disperso).

Nota metodológica. Al no cumplirse estrictamente la homogeneidad de varianzas, se realizó como verificación de robustez un ANOVA de Welch (que no asume varianzas iguales): \(F\) = 34.73, \(p\) = <1e-04. La conclusión de diferencias altamente significativas se mantiene, por lo que la violación del supuesto de homocedasticidad no compromete la conclusión general del ANOVA clásico, aunque sí sugiere tratar con algo más de cautela los valores p exactos de las comparaciones por pares que siguen a continuación (basadas en varianza combinada).

2.7 Comparaciones múltiples (Post-ANOVA) — Prueba LSD

Dado que el ANOVA resultó significativo, se procede con la prueba de Diferencia Mínima Significativa (LSD) de Fisher, adecuada aquí porque el factor tiene pocos niveles (3) definidos a priori y el interés es comparar todos los pares de medias.

pw1 <- pairwise.t.test(d1$Biomasa, d1$pH_cat, p.adjust.method = "none", pool.sd = TRUE)
medias1 <- tapply(d1$Biomasa, d1$pH_cat, mean)
ee1 <- tapply(d1$Biomasa, d1$pH_cat, function(x) sd(x)/sqrt(length(x)))
tab_lsd1 <- tabla_lsd(pw1, medias1, ee1)
kable(tab_lsd1, row.names = FALSE,
      caption = "Comparación de medias de Biomasa por nivel de pH (LSD, α = 0.05). Letras distintas = diferencia significativa")
Comparación de medias de Biomasa por nivel de pH (LSD, α = 0.05). Letras distintas = diferencia significativa
Grupo Media Letra EE
Alto 1605.39 b 141.39
Medio 1048.11 a 63.24
Bajo 593.02 c 42.56 
ggplot(d1, aes(x = pH_cat, y = Biomasa, fill = pH_cat)) +
  geom_boxplot(alpha = 0.85) +
  scale_fill_manual(values = paleta[c(3,5,1)]) +
  geom_text(data = tab_lsd1, aes(x = Grupo, y = max(d1$Biomasa) * 1.05, label = Letra),
            inherit.aes = FALSE, size = 5, fontface = "bold") +
  labs(title = "Comparaciones múltiples LSD: Biomasa según nivel de pH",
       subtitle = "Letras distintas indican diferencia significativa (α = 0.05)",
       x = "Nivel de pH", y = "Biomasa (g)") +
  theme(legend.position = "none")

Las tres letras obtenidas (a, b, c) son todas distintas, lo que indica que los tres niveles de pH difieren significativamente entre sí: el nivel Alto (media = 1605.4 g) produce significativamente más biomasa que el nivel Medio (1048.1 g), y este a su vez más que el nivel Bajo (593 g). Es decir, la relación entre pH y biomasa no solo es significativa de forma global, sino que se expresa de manera gradual y consistente en los tres tramos de acidez del suelo.

Conclusiones del Punto 1. La biomasa y sus covariables presentan comportamientos univariados coherentes con procesos biológicos/ambientales, sin anomalías destacables. De las covariables consideradas, el pH es la que más se relaciona con la biomasa (\(r \approx 0.93\)), muy por encima de la Salinidad (prácticamente sin relación) y superando también al Zinc. Al categorizar el pH en terciles, el ANOVA muestra diferencias altamente significativas en la biomasa entre niveles (\(p < 0.001\)), con un tamaño de efecto grande (\(\eta^2 \approx\) 0.59). Los residuales son normales, pero no homocedásticos, no poseen el mismo nivel de dispersión; un ANOVA de Welch confirma la misma conclusión, por lo que el hallazgo es robusto. La prueba LSD muestra que los tres niveles de pH difieren entre sí, con una relación positiva monotónica entre pH y producción de biomasa, por lo tanto se puede afirmar que suelos menos ácidos favorecen una mayor producción de esta especie forrajera.

3 Punto 2: Consumo de oxígeno en moluscos

3.1 Introducción

El consumo de oxígeno es una variable ampliamente utilizada para evaluar la actividad metabólica de los organismos, ya que el oxígeno participa en la respiración aeróbica y permite la producción de energía necesaria para procesos como el crecimiento, el movimiento, el mantenimiento celular y otras funciones fisiológicas. Por esta razón, las variaciones en el consumo de O₂ pueden utilizarse como una aproximación a los cambios en las demandas energéticas de un organismo.

En los moluscos, la respiración se encuentra estrechamente relacionada con su fisiología y con las características propias de cada especie. Además, el metabolismo puede variar entre especies debido a diferencias en su tamaño, actividad, fisiología y estrategias de obtención y utilización de energía. Para realizar comparaciones entre individuos de diferente tamaño, el consumo de oxígeno puede expresarse en relación con el peso seco del organismo, obteniendo así una medida estandarizada de consumo de O₂ por unidad de biomasa (Bayne & Newell, 1983).

3.2 Metodología

Se evaluó el consumo de oxígeno (proporción de O₂ por unidad de peso seco) de dos tipos de molusco (A y B) sometidos a tres concentraciones de agua de mar (100%, 75% y 50%), en un diseño factorial 2×3 balanceado (8 réplicas por combinación, 48 observaciones en total). Se plantea un ANOVA de dos vías con interacción: cons_o ~ molusco * c_agua, y se aplica la prueba LSD a los efectos que resulten significativos.

load("moluscos.RData")
d2 <- as.data.frame(BD_moluscos)
d2$c_agua <- factor(d2$c_agua, levels = c(50, 75, 100),
                     labels = c("50%", "75%", "100%"))
d2$molusco <- factor(d2$molusco)
kable(head(d2, 6), caption = "Primeras filas de la base de moluscos")
Primeras filas de la base de moluscos
c_agua molusco cons_o
100% A 7.16
100% A 8.26
100% A 6.78
100% A 14.00
100% A 13.60
100% A 11.10

3.3 Exploración de datos — análisis univariado

tab_u2 <- d2 %>%
  group_by(molusco, c_agua) %>%
  summarise(n = n(), Media = round(mean(cons_o),2), DE = round(sd(cons_o),2), .groups="drop")
kable(tab_u2, caption = "Estadísticos descriptivos del consumo de O2 por combinación molusco × concentración")
Estadísticos descriptivos del consumo de O2 por combinación molusco × concentración
molusco c_agua n Media DE
A 50% 8 12.18 3.09
A 75% 8 7.89 2.74
A 100% 8 9.94 2.75
B 50% 8 12.33 3.52
B 75% 8 6.10 2.74
B 100% 8 7.41 2.84 
p2a <- ggplot(d2, aes(x = cons_o, fill = molusco)) +
  geom_density(alpha = 0.55) +
  scale_fill_manual(values = paleta[c(1,2)]) +
  labs(title = "Distribución del consumo de O2 por tipo de molusco", x = "Consumo de O2", y = "Densidad")

p2b <- ggplot(d2, aes(x = molusco, y = cons_o, fill = molusco)) +
  geom_boxplot(alpha = 0.85) +
  scale_fill_manual(values = paleta[c(1,2)]) +
  labs(title = "Consumo de O2 por tipo de molusco", x = NULL, y = "Consumo de O2") +
  theme(legend.position = "none")

gridExtra::grid.arrange(p2a, p2b, ncol = 2)

La distribución del consumo de oxígeno es razonablemente simétrica y sin valores atípicos marcados, tanto para el molusco A como para el B. Comparando de forma marginal (sin considerar aún la concentración), ambos tipos de molusco presentan niveles de consumo de oxígeno similares, con medias globales de 10 (A) y 8.61 (B) — una diferencia pequeña que habrá que contrastar formalmente con el ANOVA.

3.4 Exploración de datos — análisis bivariado

ggplot(d2, aes(x = c_agua, y = cons_o, fill = molusco)) +
  geom_boxplot(alpha = 0.85, position = position_dodge(0.75)) +
  scale_fill_manual(values = paleta[c(1,2)]) +
  labs(title = "Consumo de O2 según concentración de agua de mar y tipo de molusco",
       x = "Concentración de agua de mar", y = "Consumo de O2", fill = "Molusco")

d2_medias <- d2 %>% group_by(molusco, c_agua) %>% summarise(media = mean(cons_o), ee = sd(cons_o)/sqrt(n()), .groups="drop")

ggplot(d2_medias, aes(x = c_agua, y = media, color = molusco, group = molusco)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = media - ee, ymax = media + ee), width = 0.1) +
  scale_color_manual(values = paleta[c(1,2)]) +
  labs(title = "Gráfico de interacción: molusco × concentración de agua de mar",
       x = "Concentración de agua de mar", y = "Consumo medio de O2 (± EE)", color = "Molusco")

En ambos tipos de molusco, el patrón general es el mismo: el consumo de oxígeno es máximo al 50% de agua de mar y disminuye al aumentar la concentración (75% y 100%), aunque la caída entre 75% y 100% es más leve e incluso se invierte levemente en el molusco A. Las líneas del gráfico de interacción son casi paralelas, lo que sugiere, de forma exploratoria, que el efecto de la concentración es similar para los dos tipos de molusco, es decir, no hay evidencia visual fuerte de interacción — hipótesis que se confirma a continuación con el ANOVA.

3.5 Análisis de Varianza de dos vías

\[cons\_o_{ijk} = \mu + \alpha_i(\text{molusco}) + \beta_j(\text{concentración}) + (\alpha\beta)_{ij} + \varepsilon_{ijk}\]

mod2 <- aov(cons_o ~ molusco * c_agua, data = d2)
tabla_anova2 <- summary(mod2)[[1]]
kable(round(tabla_anova2, 4), caption = "Tabla ANOVA de dos vías — Consumo de O2 ~ Molusco × Concentración")
Tabla ANOVA de dos vías — Consumo de O2 ~ Molusco × Concentración
Df Sum Sq Mean Sq F value Pr(>F)
molusco 1 23.2269 23.2269 2.6508 0.1110
c_agua 2 230.8160 115.4080 13.1712 0.0000
molusco:c_agua 2 15.3563 7.6781 0.8763 0.4238
Residuals 42 368.0112 8.7622 NA NA

Resultados.

  • Molusco (efecto principal): \(F\)(NA,42) = NA, \(p\) = NA → no significativo. No hay evidencia de que el tipo de molusco, por sí solo, afecte el consumo de oxígeno.
  • Concentración de agua de mar (efecto principal): \(F\)(2,42) = 13.17, \(p\) = <1e-04 → altamente significativo. La concentración de agua de mar sí afecta el consumo de oxígeno.
  • Interacción molusco:concentración: \(F\)(2,42) = 0.88, \(p\) = 0.424 → no significativa. El efecto de la concentración es estadísticamente el mismo para ambos tipos de molusco (consistente con las líneas casi paralelas del gráfico de interacción).

3.6 Verificación de supuestos

p2q <- ggplot(data.frame(res = residuals(mod2)), aes(sample = res)) +
  stat_qq(color = paleta[1]) + stat_qq_line(color = "#B3703A") +
  labs(title = "Q-Q plot de los residuales", x="Cuantiles teóricos", y="Cuantiles muestrales")
p2r <- ggplot(data.frame(res = residuals(mod2), fit = fitted(mod2)), aes(x = fit, y = res)) +
  geom_point(color = paleta[1]) + geom_hline(yintercept = 0, linetype="dashed", color="#B3703A") +
  labs(title = "Residuales vs. valores ajustados", x = "Valores ajustados", y = "Residuales")
gridExtra::grid.arrange(p2q, p2r, ncol = 2)

sh2 <- shapiro.test(residuals(mod2))
lv2 <- leveneTest(cons_o ~ molusco * c_agua, data = d2)
kable(data.frame(Prueba = c("Shapiro-Wilk (normalidad residuales)", "Levene (homogeneidad de varianzas)"),
                  Estadístico = c(round(sh2$statistic,3), round(lv2$`F value`[1],3)),
                  Valor_p = c(format.pval(sh2$p.value, digits=3),
                              format.pval(lv2$`Pr(>F)`[1], digits=3))),
      caption = "Verificación de supuestos del modelo ANOVA — Punto 2")
Verificación de supuestos del modelo ANOVA — Punto 2
Prueba Estadístico Valor_p
W Shapiro-Wilk (normalidad residuales) 0.958 0.0857
Levene (homogeneidad de varianzas) 0.172 0.971

Tanto la prueba de Shapiro-Wilk (\(p\) = 0.0857) como la de Levene (\(p\) = 0.971) arrojan valores \(p\) superiores a 0.05, por lo que no se rechazan los supuestos de normalidad de los residuales ni de homogeneidad de varianzas. El modelo ANOVA de dos vías es, por tanto, válido para este conjunto de datos, y sus conclusiones pueden interpretarse con confianza.

3.7 Comparaciones múltiples (Post-ANOVA) — Prueba LSD

Como el factor molusco y la interacción no fueron significativos, la comparación post-hoc solo tiene sentido para el factor concentración de agua de mar, que sí resultó significativo.

pw2 <- pairwise.t.test(d2$cons_o, d2$c_agua, p.adjust.method = "none", pool.sd = TRUE)
medias2 <- tapply(d2$cons_o, d2$c_agua, mean)
ee2 <- tapply(d2$cons_o, d2$c_agua, function(x) sd(x)/sqrt(length(x)))
tab_lsd2 <- tabla_lsd(pw2, medias2, ee2)
kable(tab_lsd2, row.names = FALSE,
      caption = "Comparación de medias de consumo de O2 por concentración de agua de mar (LSD, α = 0.05)")
Comparación de medias de consumo de O2 por concentración de agua de mar (LSD, α = 0.05)
Grupo Media Letra EE
50% 12.25 b 0.80
100% 8.67 a 0.75
75% 6.99 a 0.70 
ggplot(d2, aes(x = c_agua, y = cons_o, fill = c_agua)) +
  geom_boxplot(alpha = 0.85) +
  scale_fill_manual(values = paleta[c(3,2,1)]) +
  geom_text(data = tab_lsd2, aes(x = Grupo, y = max(d2$cons_o) * 1.05, label = Letra),
            inherit.aes = FALSE, size = 5, fontface = "bold") +
  labs(title = "Comparaciones múltiples LSD: consumo de O2 por concentración",
       subtitle = "Letras distintas indican diferencia significativa (α = 0.05)",
       x = "Concentración de agua de mar", y = "Consumo de O2") +
  theme(legend.position = "none")

La concentración de 50% (media = 12.25) presenta un consumo de oxígeno significativamente mayor que las concentraciones de 75% (6.99, letra distinta) y de 100% (8.67, letra distinta). Sin embargo, 75% y 100% no difieren significativamente entre sí (comparten letra), es decir, una vez que la concentración de agua de mar baja de 100%, el mayor incremento en el consumo de oxígeno ocurre al llegar al 50%, mientras que la diferencia entre concentraciones intermedias y altas no alcanza a ser estadísticamente relevante.

Conclusiones del Punto 2. El consumo de oxígeno no difiere significativamente entre los dos tipos de molusco (A y B) considerados de forma aislada. La concentración de agua de mar sí es un factor determinante del consumo de oxígeno (\(p<0.001\)). No existe interacción significativa entre tipo de molusco y concentración: el patrón de respuesta a la concentración es el mismo para ambas especies (líneas prácticamente paralelas). Los supuestos de normalidad y homocedasticidad se cumplen, por lo que el modelo es confiable. La prueba LSD indica que el consumo de oxígeno se dispara ante condiciones de estrés osmótico (agua de mar diluida al 50%), mientras que concentraciones de 75% y 100% generan un metabolismo respiratorio comparable entre sí.

4 Punto 3: Biodiversidad de anfibios y uso del suelo

4.1 Introducción

Los anfibios constituyen un grupo de vertebrados estrechamente relacionado con las condiciones ambientales de los ecosistemas debido a características de su ciclo de vida y de su fisiología. Su dependencia de ambientes húmedos y su sensibilidad a modificaciones en las condiciones del hábitat hacen que sean considerados un grupo de interés para evaluar cambios en los ecosistemas y posibles efectos de la intervención humana sobre la biodiversidad.

La transformación de los bosques naturales para actividades productivas puede modificar la estructura de la vegetación, la disponibilidad de refugios y sitios de reproducción y las condiciones microclimáticas, generando cambios en la composición de las comunidades de anfibios. Por ello, comparar ambientes con diferentes niveles de intervención permite evaluar cómo la modificación del uso del suelo se relaciona con la biodiversidad presente en una reserva forestal.

Adicionalmente, el concepto de riqueza de especies corresponde al número de especies diferentes presentes en una comunidad o área determinada, por lo que constituye una medida básica de la diversidad biológica. Sin embargo, esta medida no considera cuántos individuos pertenecen a cada especie. Para incorporar tanto la riqueza como la distribución de las abundancias entre las especies, se utilizan índices de diversidad, entre ellos el índice de Shannon-Wiener (H’). Este índice permite caracterizar la diversidad de una comunidad teniendo en cuenta simultáneamente el número de especies presentes y qué tan uniformemente se distribuyen sus individuos entre ellas (Magurran,2004). En consecuencia, comunidades con mayor riqueza y una distribución más equitativa de las abundancias tienden a presentar valores más altos de diversidad de Shannon-Wiener.

4.2 Metodología

Se evaluó el efecto del tipo de hábitat (cuatro niveles que representan un gradiente de intervención antrópica: Bosque primario, Bosque secundario, Sistema silvopastoril y Potrero) sobre la riqueza de especies, el índice de diversidad de Shannon-Wiener y se registró además la altitud de cada una de las 52 parcelas (13 por hábitat). El literal (c) pide evaluar mediante un ANOVA de una vía si el hábitat afecta significativamente el índice de Shannon, con LSD como prueba post-hoc.

load("Biodiversidad.RData")
d3 <- BD_biodiversidad
d3$Habitat <- factor(d3$Habitat,
                      levels = c("Bosque primario","Bosque secundario","Sistema silvopastoril","Potrero"))
kable(head(d3, 6), caption = "Primeras filas de la base de biodiversidad")
Primeras filas de la base de biodiversidad
Parcela Habitat Riqueza Shannon Altitud
P001 Bosque primario 14 2.13 1478.0
P002 Bosque primario 14 2.98 1476.8
P003 Bosque primario 13 2.07 1505.2
P004 Bosque primario 14 2.53 1395.2
P005 Bosque primario 10 2.82 1493.5
P006 Bosque primario 15 2.75 1337.9

4.3 Exploración de datos — análisis univariado

tab_u3 <- rbind(
  resumen_num(d3$Riqueza, "Riqueza (n° especies)"),
  resumen_num(d3$Shannon, "Índice de Shannon"),
  resumen_num(d3$Altitud, "Altitud (m.s.n.m.)")
)
kable(tab_u3, caption = "Estadísticos descriptivos generales (n = 52 parcelas)")
Estadísticos descriptivos generales (n = 52 parcelas)
Variable n Media DE CV_pct Mínimo Q1 Mediana Q3 Máximo Asimetría
25% Riqueza (n° especies) 52 9.73 4.37 44.9 2.00 6.00 9.50 13.00 21.00 0.28
25%1 Índice de Shannon 52 1.76 0.61 34.4 0.67 1.26 1.85 2.11 2.98 0.04
25%2 Altitud (m.s.n.m.) 52 1243.89 172.71 13.9 950.90 1081.88 1254.80 1382.88 1600.80 0.09 
wrap_plots(
  poligono_box(d3$Riqueza, "Riqueza de especies", paleta[1]),
  poligono_box(d3$Shannon, "Índice de Shannon",   paleta[2]),
  poligono_box(d3$Altitud, "Altitud (m.s.n.m.)",  paleta[4]),
  ncol = 3
) + plot_annotation(title = "Distribución univariada de riqueza, diversidad y altitud",
                     theme = theme(plot.title = element_text(face="bold", color="#3B5D42", family="serif")))

La riqueza de especies varía entre 2 y 21 especies por parcela, con una distribución razonablemente simétrica. El índice de Shannon se mueve entre 0.67 y 2.98, valores consistentes con comunidades de baja a moderada diversidad. La altitud presenta un rango amplio (951–1601 m.s.n.m.), reflejo de que las parcelas de los distintos hábitats fueron establecidas en zonas topográficamente diferentes de la reserva.

4.4 Exploración de datos — análisis bivariado

p3d <- ggplot(d3, aes(x = Habitat, y = Riqueza, fill = Habitat)) +
  geom_boxplot(alpha = 0.85) +
  scale_fill_manual(values = paleta) +
  labs(title = "Riqueza por tipo de hábitat", x = NULL, y = "Riqueza") +
  theme(legend.position = "none", axis.text.x = element_text(angle = 25, hjust = 1))

p3e <- ggplot(d3, aes(x = Habitat, y = Shannon, fill = Habitat)) +
  geom_boxplot(alpha = 0.85) +
  scale_fill_manual(values = paleta) +
  labs(title = "Shannon por tipo de hábitat", x = NULL, y = "Índice de Shannon") +
  theme(legend.position = "none", axis.text.x = element_text(angle = 25, hjust = 1))

gridExtra::grid.arrange(p3d, p3e, ncol = 2)

cor_alt <- cor(d3$Altitud, d3$Shannon)
ggplot(d3, aes(x = Altitud, y = Shannon, color = Habitat)) +
  geom_point(size = 2.4, alpha = 0.85) +
  geom_smooth(aes(group = 1), method = "lm", se = TRUE, color = "#3B5D42", linetype = "dashed") +
  scale_color_manual(values = paleta) +
  labs(title = "Relación entre altitud y diversidad de Shannon",
       subtitle = paste0("Correlación de Pearson global: r = ", round(cor_alt, 2)),
       x = "Altitud (m.s.n.m.)", y = "Índice de Shannon")

Tanto la riqueza como el índice de Shannon muestran un patrón decreciente y muy marcado a lo largo del gradiente de intervención: Bosque primario > Bosque secundario > Sistema silvopastoril > Potrero, con muy poco solapamiento entre las cajas de los hábitats extremos. En cuanto a la relación con la altitud, se observa una correlación positiva fuerte (\(r =\) 0.75) entre altitud y diversidad de Shannon; sin embargo, esto se explica en buena parte porque los hábitats menos intervenidos (bosques) de esta reserva se ubican también en las zonas más altas, por lo que la altitud actúa aquí como una variable asociada al hábitat más que como una causa independiente de la diversidad (Busch & Ferretti-Gallon, 2023).

4.5 Análisis de Varianza de una vía

\[Shannon_{ij} = \mu + \tau_i(\text{Hábitat}) + \varepsilon_{ij}, \qquad i = 1,\dots,4\]

mod3 <- aov(Shannon ~ Habitat, data = d3)
tabla_anova3 <- summary(mod3)[[1]]
kable(round(tabla_anova3, 4), caption = "Tabla ANOVA — Índice de Shannon según tipo de hábitat")
Tabla ANOVA — Índice de Shannon según tipo de hábitat
Df Sum Sq Mean Sq F value Pr(>F)
Habitat 3 14.8624 4.9541 60.6112 0
Residuals 48 3.9233 0.0817 NA NA 
eta2_3 <- tabla_anova3["Habitat","Sum Sq"] / sum(tabla_anova3[,"Sum Sq"])

Resultado del ANOVA. \(F\)(3,48) = 60.61, \(p\) = <1e-04. Se rechaza contundentemente la hipótesis nula: el tipo de hábitat sí genera diferencias muy significativas en el índice de diversidad de Shannon. El tamaño del efecto es muy grande (\(\eta^2 =\) 0.79): el hábitat explica cerca del 79% de la variación total observada en la diversidad.

4.6 Verificación de supuestos

p3q <- ggplot(data.frame(res = residuals(mod3)), aes(sample = res)) +
  stat_qq(color = paleta[1]) + stat_qq_line(color = "#B3703A") +
  labs(title = "Q-Q plot de los residuales", x="Cuantiles teóricos", y="Cuantiles muestrales")
p3r <- ggplot(data.frame(res = residuals(mod3), fit = fitted(mod3)), aes(x = fit, y = res)) +
  geom_point(color = paleta[1]) + geom_hline(yintercept = 0, linetype="dashed", color="#B3703A") +
  labs(title = "Residuales vs. valores ajustados", x = "Valores ajustados", y = "Residuales")
gridExtra::grid.arrange(p3q, p3r, ncol = 2)

sh3 <- shapiro.test(residuals(mod3))
lv3 <- leveneTest(Shannon ~ Habitat, data = d3)
kable(data.frame(Prueba = c("Shapiro-Wilk (normalidad residuales)", "Levene (homogeneidad de varianzas)"),
                  Estadístico = c(round(sh3$statistic,3), round(lv3$`F value`[1],3)),
                  Valor_p = c(format.pval(sh3$p.value, digits=3),
                              format.pval(lv3$`Pr(>F)`[1], digits=3))),
      caption = "Verificación de supuestos del modelo ANOVA — Punto 3")
Verificación de supuestos del modelo ANOVA — Punto 3
Prueba Estadístico Valor_p
W Shapiro-Wilk (normalidad residuales) 0.978 0.463
Levene (homogeneidad de varianzas) 1.277 0.293

Los valores \(p\) de Shapiro-Wilk (\(p\) = 0.463) y de Levene (\(p\) = 0.293) son ambos superiores a 0.05: no hay evidencia para rechazar ni la normalidad de los residuales ni la homogeneidad de varianzas entre hábitats. El modelo ANOVA de una vía es, por tanto, plenamente válido para estos datos.

4.7 Comparaciones múltiples (Post-ANOVA) — Prueba LSD

pw3 <- pairwise.t.test(d3$Shannon, d3$Habitat, p.adjust.method = "none", pool.sd = TRUE)
medias3 <- tapply(d3$Shannon, d3$Habitat, mean)
ee3 <- tapply(d3$Shannon, d3$Habitat, function(x) sd(x)/sqrt(length(x)))
tab_lsd3 <- tabla_lsd(pw3, medias3, ee3)
kable(tab_lsd3, row.names = FALSE,
      caption = "Comparación de medias del índice de Shannon por hábitat (LSD, α = 0.05)")
Comparación de medias del índice de Shannon por hábitat (LSD, α = 0.05)
Grupo Media Letra EE
Bosque primario 2.47 d 0.10
Bosque secundario 1.97 a 0.06
Sistema silvopastoril 1.60 b 0.08
Potrero 1.00 c 0.07 
ggplot(d3, aes(x = Habitat, y = Shannon, fill = Habitat)) +
  geom_boxplot(alpha = 0.85) +
  scale_fill_manual(values = paleta) +
  geom_text(data = tab_lsd3, aes(x = Grupo, y = max(d3$Shannon) * 1.08, label = Letra),
            inherit.aes = FALSE, size = 5, fontface = "bold") +
  labs(title = "Comparaciones múltiples LSD: índice de Shannon por hábitat",
       subtitle = "Letras distintas indican diferencia significativa (α = 0.05)",
       x = NULL, y = "Índice de Shannon") +
  theme(legend.position = "none", axis.text.x = element_text(angle = 20, hjust = 1))

Las cuatro letras obtenidas son todas diferentes entre sí (a, b, c, d), es decir, cada tipo de hábitat difiere significativamente de todos los demás en su índice de diversidad de Shannon. El orden de las medias reproduce exactamente el gradiente de intervención antrópica: Bosque primario (2.47) > Bosque secundario (1.97) > Sistema silvopastoril (1.6) > Potrero (1). No existen, por tanto, hábitats “empatados” en diversidad: cada nivel de intervención adicional sobre el bosque original conlleva una pérdida estadísticamente detectable de diversidad de anfibios.

Conclusiones del Punto 3. La riqueza y la diversidad de anfibios disminuyen de forma marcada y consistente a lo largo del gradiente de intervención antrópica del suelo. El ANOVA de una vía confirma diferencias muy significativas en el índice de Shannon entre hábitats (\(p<0.001\)), con un tamaño de efecto muy grande (\(\eta^2 \approx\) 0.79). Los supuestos de normalidad y homocedasticidad se cumplen sin reservas, lo que da robustez estadística a la conclusión. La prueba LSD muestra que los cuatro hábitats son distinguibles entre sí: el paso de bosque primario a potrero implica una pérdida progresiva y significativa de diversidad en cada etapa del gradiente. La fuerte asociación entre altitud y diversidad observada en el análisis bivariado debe interpretarse con cautela, pues está confundida con el tipo de hábitat (los bosques mejor conservados de esta reserva se ubican en las zonas más altas).

5 Conclusiones generales del taller

A través de los tres estudios de caso se evidencia el valor del ANOVA como herramienta para formalizar diferencias que, en la exploración descriptiva, ya se insinuaban gráficamente:

  • En el Punto 1, el análisis exploratorio bivariado fue clave para identificar correctamente la covariable relevante (pH, y no la Salinidad que da nombre a la base) antes de construir el factor categórico del ANOVA.
  • En el Punto 2, el diseño factorial permitió separar el efecto de la concentración de agua de mar (relevante) del efecto del tipo de molusco (no relevante) y descartar formalmente una interacción que parecía poco probable visualmente.
  • En el Punto 3, el ANOVA de una vía confirmó de manera contundente que el uso del suelo es un determinante fuerte de la biodiversidad de anfibios, con diferencias significativas entre todos los pares de hábitats.

6 Referencias

  • Montgomery, D. C. (2017). Design and Analysis of Experiments (9th ed.). Wiley.
  • Fox, J. & Weisberg, S. (2019). An R Companion to Applied Regression (3rd ed.). Sage.
  • Neina, D. (2019). The Role of Soil pH in Plant Nutrition and Soil Remediation. Applied and Environmental Soil Science, 2019, Article 5794869.
  • Bayne, B. L., & Newell, R. C. (1983). Physiological energetics of marine molluscs. En A. S. M. Saleuddin & K. M. Wilbur (Eds.), The Mollusca: Physiology, Part 1 (Vol. 4, pp. 407–515). Academic Press.
  • Magurran, A. E. (2004). Measuring Biological Diversity. Blackwell Publishing.
  • Busch, J., & Ferretti-Gallon, K. (2023). What drives and stops deforestation, reforestation, and forest degradation? An updated meta-analysis. Review of Environmental Economics and Policy, 17(2).