Este informe desarrolla los tres puntos del Taller 1 de Diseño de Experimentos: análisis exploratorio, ANOVA y pruebas post-ANOVA (LSD de Fisher) para tres conjuntos de datos independientes: Salinidad (biomasa vegetal y covariables de suelo), Moluscos (consumo de oxígeno bajo distintas concentraciones de agua de mar) y Biodiversidad (riqueza y diversidad de anfibios según el uso del suelo).
# Carga reproducible de las tres bases de datos
load("Salinidad.RData") # objeto: Salinidad
load("moluscos.RData") # objeto: BD_moluscos
load("Biodiversidad.RData") # objeto: BD_biodiversidad
BD_moluscos <- as.data.frame(BD_moluscos)
BD_moluscos$molusco <- factor(BD_moluscos$molusco)
BD_moluscos$c_agua <- factor(BD_moluscos$c_agua, levels = c(50, 75, 100))
BD_biodiversidad$Habitat <- factor(BD_biodiversidad$Habitat,
levels = c("Bosque primario", "Bosque secundario", "Sistema silvopastoril", "Potrero"))resumen1 <- Salinidad %>%
summarise(across(everything(), list(Media = mean, DE = sd, Min = min, Max = max, CV = ~sd(.)/mean(.)*100))) %>%
pivot_longer(everything(), names_to = "var", values_to = "valor") %>%
separate(var, into = c("Variable", "Estadístico"), sep = "_(?=[^_]+$)") %>%
pivot_wider(names_from = Estadístico, values_from = valor)
tabla_bonita(resumen1, digits = 2, caption = "Estadísticos descriptivos — datos Salinidad")| Variable | Media | DE | Min | Max | CV |
|---|---|---|---|---|---|
| Biomasa | 1082.17 | 546.29 | 369.82 | 2337.33 | 50.48 |
| pH | 4.61 | 1.25 | 3.20 | 7.45 | 27.22 |
| Salinidad | 30.27 | 3.72 | 24.00 | 38.00 | 12.29 |
| Zinc | 17.83 | 8.27 | 0.21 | 31.29 | 46.40 |
| Potasio | 797.38 | 297.58 | 350.73 | 1441.67 | 37.32 |
vars1 <- c("Biomasa", "pH", "Salinidad", "Zinc", "Potasio")
histos <- lapply(vars1, function(v) {
ggplot(Salinidad, aes(x = .data[[v]])) +
geom_histogram(fill = morado, color = "white", bins = 10, alpha = 0.85) +
labs(title = v, x = v, y = "Frecuencia") +
theme_michi(10)
})
do.call(grid.arrange, c(histos, ncol = 3))Interpretación: la biomasa presenta el mayor coeficiente de variación entre las variables medidas, reflejando la heterogeneidad ambiental de las 45 muestras. El pH se concentra entre 4 y 5.5 (rango típico de suelos ácidos), la salinidad varía poco (CV bajo), y tanto el zinc como el potasio muestran dispersión moderada-alta, sin evidencia fuerte de asimetría extrema en ninguna variable según los histogramas.
p_ph <- ggplot(Salinidad, aes(pH, Biomasa)) + geom_point(color = azul, size = 2) +
geom_smooth(method = "lm", color = morado, se = FALSE) + theme_michi(10) + labs(title = "Biomasa vs pH")
p_sal <- ggplot(Salinidad, aes(Salinidad, Biomasa)) + geom_point(color = azul, size = 2) +
geom_smooth(method = "lm", color = morado, se = FALSE) + theme_michi(10) + labs(title = "Biomasa vs Salinidad")
p_zn <- ggplot(Salinidad, aes(Zinc, Biomasa)) + geom_point(color = azul, size = 2) +
geom_smooth(method = "lm", color = morado, se = FALSE) + theme_michi(10) + labs(title = "Biomasa vs Zinc")
grid.arrange(p_ph, p_sal, p_zn, ncol = 3)correlaciones <- cor(Salinidad$Biomasa, Salinidad[, c("pH", "Salinidad", "Zinc", "Potasio")])
tabla_bonita(round(correlaciones, 3), caption = "Correlación de Pearson con Biomasa")| pH | Salinidad | Zinc | Potasio |
|---|---|---|---|
| 0.928 | -0.067 | -0.781 | -0.073 |
Interpretación: el pH presenta la correlación más fuerte con la biomasa (r ≈ 0.93, relación positiva casi lineal), seguido por el Zinc (r ≈ −0.78, relación negativa moderada-fuerte). La salinidad y el potasio muestran correlaciones prácticamente nulas con la biomasa en este conjunto de datos. Por lo tanto, el pH es la covariable con mayor relación con la biomasa y será la variable categorizada en el literal (c).
terciles_ph <- quantile(Salinidad$pH, probs = c(1/3, 2/3))
Salinidad$pH_cat <- cut(Salinidad$pH,
breaks = c(-Inf, terciles_ph[1], terciles_ph[2], Inf),
labels = c("Bajo", "Medio", "Alto"))
kable(table(Salinidad$pH_cat), col.names = c("Nivel de pH", "n")) %>%
kable_styling(full_width = FALSE)| Nivel de pH | n |
|---|---|
| Bajo | 15 |
| Medio | 15 |
| Alto | 15 |
Se usaron terciles (percentiles 33 % y 67 %) como criterio de corte porque generan grupos de tamaño balanceado (15 observaciones cada uno), lo cual favorece la validez del ANOVA y de la prueba LSD posterior.
ggplot(Salinidad, aes(pH_cat, Biomasa, fill = pH_cat)) +
geom_boxplot(alpha = 0.85) +
scale_fill_michi() +
labs(title = "Biomasa según nivel de pH del suelo", x = "Nivel de pH", y = "Biomasa (g)", fill = "Nivel") +
theme_michi()mod1 <- aov(Biomasa ~ pH_cat, data = Salinidad)
tabla_bonita(as.data.frame(summary(mod1)[[1]]), digits = 3, caption = "Tabla ANOVA — Biomasa ~ Nivel de pH")| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| pH_cat | 2 | 7712683 | 3856341.6 | 29.893 | 0 |
| Residuals | 42 | 5418235 | 129005.6 | NA | NA |
res1 <- residuals(mod1)
sw1 <- shapiro.test(res1)
lev1 <- leveneTest(Biomasa ~ pH_cat, data = Salinidad)
p_qq <- ggplot(data.frame(res1), aes(sample = res1)) +
stat_qq(color = azul) + stat_qq_line(color = morado) +
labs(title = "QQ-plot de residuales") + theme_michi(10)
p_res <- ggplot(data.frame(fit = fitted(mod1), res = res1), aes(fit, res)) +
geom_point(color = azul) + geom_hline(yintercept = 0, color = morado, linetype = "dashed") +
labs(title = "Residuales vs ajustados", x = "Ajustados", y = "Residuales") + theme_michi(10)
grid.arrange(p_qq, p_res, ncol = 2)El supuesto de normalidad se cumple, pero el de homocedasticidad se viola (Levene p < 0.001) — el grupo “Alto” tiene una desviación estándar (≈548) muy superior a la del grupo “Bajo” (≈165). Como verificación de robustez ante esta violación, se ejecuta una prueba de Welch (no asume varianzas iguales):
##
## One-way analysis of means (not assuming equal variances)
##
## data: Biomasa and pH_cat
## F = 34.727, num df = 2.000, denom df = 24.736, p-value = 6.58e-08
La conclusión de significancia se mantiene bajo Welch (p < 0.001), por lo que se procede con la prueba post-hoc, interpretando la LSD con cautela dado el supuesto violado.
lsd1 <- lsd_test(Salinidad$Biomasa, Salinidad$pH_cat, mod1)
tabla_bonita(lsd1$comparaciones, caption = "Comparaciones LSD por pares")| Comparacion | Diferencia | Valor_p | Significativo | |
|---|---|---|---|---|
| Bajo | Bajo vs Medio | -455.086 | 0.0012 | Si |
| Bajo1 | Bajo vs Alto | -1012.362 | 0.0000 | Si |
| Medio | Medio vs Alto | -557.276 | 0.0001 | Si |
tabla_bonita(lsd1$resumen, caption = "Medias y grupos de significancia (letras distintas = diferencia significativa)")| Grupo | Media | n | Letra | |
|---|---|---|---|---|
| Alto | Alto | 1605.385 | 15 | c |
| Medio | Medio | 1048.109 | 15 | b |
| Bajo | Bajo | 593.023 | 15 | a |
Interpretación: el ANOVA detecta diferencias altamente significativas en la biomasa según el nivel de pH (p < 0.001). La prueba LSD muestra que los tres niveles difieren significativamente entre sí: a mayor pH del suelo, mayor biomasa producida (Bajo ≈ 593 g, Medio ≈ 1048 g, Alto ≈ 1605 g), consistente con la fuerte correlación positiva observada en el literal (b). Suelos menos ácidos favorecen notoriamente la producción de biomasa de esta especie forrajera.
resumen2 <- BD_moluscos %>%
group_by(molusco, c_agua) %>%
summarise(Media = mean(cons_o), DE = sd(cons_o), Min = min(cons_o), Max = max(cons_o), n = n(), .groups = "drop")
tabla_bonita(resumen2, digits = 2, caption = "Consumo de O2 por tipo de molusco y concentración de agua de mar")| molusco | c_agua | Media | DE | Min | Max | n |
|---|---|---|---|---|---|---|
| A | 50 | 12.18 | 3.09 | 9.74 | 18.80 | 8 |
| A | 75 | 7.89 | 2.74 | 5.20 | 13.20 | 8 |
| A | 100 | 9.94 | 2.75 | 6.78 | 14.00 | 8 |
| B | 50 | 12.33 | 3.52 | 6.38 | 17.70 | 8 |
| B | 75 | 6.10 | 2.74 | 1.80 | 9.96 | 8 |
| B | 100 | 7.41 | 2.84 | 3.68 | 11.60 | 8 |
h1 <- ggplot(BD_moluscos, aes(cons_o)) + geom_histogram(fill = cian, color = "white", bins = 10) +
labs(title = "Distribución del consumo de O2") + theme_michi(10)
h2 <- ggplot(BD_moluscos, aes(molusco, cons_o, fill = molusco)) + geom_boxplot() +
scale_fill_michi() + labs(title = "Consumo de O2 por tipo de molusco") + theme_michi(10)
grid.arrange(h1, h2, ncol = 2)Interpretación: el consumo de oxígeno se distribuye de forma razonablemente simétrica, sin valores atípicos marcados. Los dos tipos de molusco muestran rangos y medianas similares vistos de forma marginal, lo que sugiere que la diferencia entre ellos (si existe) podría depender de la concentración de agua de mar.
b1 <- ggplot(BD_moluscos, aes(c_agua, cons_o, fill = c_agua)) + geom_boxplot(alpha = 0.85) +
scale_fill_michi() + labs(title = "Consumo de O2 según concentración de agua de mar",
x = "Concentración de agua de mar (%)", y = "Consumo de O2") + theme_michi(10)
resumen_int <- BD_moluscos %>% group_by(molusco, c_agua) %>% summarise(media = mean(cons_o), .groups = "drop")
b2 <- ggplot(resumen_int, aes(c_agua, media, color = molusco, group = molusco)) +
geom_point(size = 3) + geom_line(linewidth = 1) +
scale_color_michi() +
labs(title = "Interacción molusco × concentración", x = "Concentración (%)", y = "Consumo medio de O2", color = "Molusco") +
theme_michi(10)
grid.arrange(b1, b2, ncol = 2)Interpretación: el consumo de oxígeno disminuye claramente al pasar de 50 % a 75-100 % de concentración de agua de mar, en ambos tipos de molusco. El gráfico de interacción muestra líneas aproximadamente paralelas entre A y B, lo que sugiere que el patrón de respuesta a la concentración de agua de mar es similar para ambos tipos de molusco (sin evidencia visual fuerte de interacción). Sin embargo, el patrón no es estrictamente decreciente: el consumo pasa de 12.25 mgO₂ (50 %) a 6.99 mgO₂ (75 %) y luego sube ligeramente a 8.67 mgO₂ (100 %), aunque esta diferencia entre 75 % y 100 % no resultó significativa (LSD, p = 0.12). Esto es consistente con la fisiología osmorreguladora de los moluscos marinos: el agua de mar al 100 % es prácticamente isosmótica con los fluidos internos del organismo, por lo que no exige gasto energético adicional, mientras que diluir el agua de mar somete al molusco a estrés osmótico y obliga a invertir energía en osmorregulación, lo que eleva su tasa metabólica y, por tanto, su consumo de oxígeno. El leve repunte no significativo en 75 % sugiere que la relación entre dilución y gasto osmorregulador no es estrictamente lineal, posiblemente por un umbral de tolerancia antes de que el estrés se vuelva metabólicamente costoso.
mod2 <- aov(cons_o ~ molusco * c_agua, data = BD_moluscos)
tabla_bonita(as.data.frame(summary(mod2)[[1]]), digits = 4, caption = "Tabla ANOVA de dos vías — Consumo de O2")| 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 |
res2 <- residuals(mod2)
sw2 <- shapiro.test(res2)
lev2 <- leveneTest(cons_o ~ molusco*c_agua, data = BD_moluscos)
p_qq2 <- ggplot(data.frame(res2), aes(sample = res2)) + stat_qq(color = azul) + stat_qq_line(color = morado) +
labs(title = "QQ-plot de residuales") + theme_michi(10)
p_res2 <- ggplot(data.frame(fit = fitted(mod2), res = res2), aes(fit, res)) +
geom_point(color = azul) + geom_hline(yintercept = 0, color = morado, linetype = "dashed") +
labs(title = "Residuales vs ajustados") + theme_michi(10)
grid.arrange(p_qq2, p_res2, ncol = 2)Ambos supuestos se satisfacen, por lo que el ANOVA de dos vías y la prueba LSD son directamente interpretables sin necesidad de alternativas robustas.
lsd2 <- lsd_test(BD_moluscos$cons_o, BD_moluscos$c_agua, mod2)
tabla_bonita(lsd2$comparaciones, caption = "Comparaciones LSD — concentración de agua de mar")| Comparacion | Diferencia | Valor_p | Significativo | |
|---|---|---|---|---|
| 50 | 50 vs 75 | 5.258 | 0.0000 | Si |
| 501 | 50 vs 100 | 3.579 | 0.0014 | Si |
| 75 | 75 vs 100 | -1.679 | 0.1162 | No |
| Grupo | Media | n | Letra | |
|---|---|---|---|---|
| 50 | 50 | 12.251 | 16 | a |
| 100 | 100 | 8.671 | 16 | b |
| 75 | 75 | 6.992 | 16 | b |
Interpretación: dado que solo el factor concentración resultó significativo, la LSD se aplicó sobre sus tres niveles. El consumo de oxígeno a 50 % es significativamente mayor que a 75 % y 100 %, mientras que 75 % y 100 % no difieren significativamente entre sí. Esto sugiere que la mayor disponibilidad de agua dulce (menor concentración de agua de mar) incrementa la actividad metabólica de ambos moluscos, posiblemente por estrés osmótico en concentraciones altas.
resumen3 <- BD_biodiversidad %>%
summarise(across(c(Riqueza, Shannon, Altitud),
list(Media = mean, DE = sd, Min = min, Max = max))) %>%
pivot_longer(everything(), names_to = "var", values_to = "valor") %>%
separate(var, into = c("Variable", "Estadístico"), sep = "_(?=[^_]+$)") %>%
pivot_wider(names_from = Estadístico, values_from = valor)
tabla_bonita(resumen3, digits = 2, caption = "Estadísticos descriptivos — riqueza, Shannon y altitud")| Variable | Media | DE | Min | Max |
|---|---|---|---|---|
| Riqueza | 9.73 | 4.37 | 2.00 | 21.00 |
| Shannon | 1.76 | 0.61 | 0.67 | 2.98 |
| Altitud | 1243.89 | 172.71 | 950.90 | 1600.80 |
h3a <- ggplot(BD_biodiversidad, aes(Riqueza)) + geom_histogram(fill = azul, color = "white", bins = 8) +
theme_michi(10) + labs(title = "Riqueza de especies")
h3b <- ggplot(BD_biodiversidad, aes(Shannon)) + geom_histogram(fill = morado, color = "white", bins = 8) +
theme_michi(10) + labs(title = "Índice de Shannon-Wiener")
h3c <- ggplot(BD_biodiversidad, aes(Altitud)) + geom_histogram(fill = cian, color = "white", bins = 8) +
theme_michi(10) + labs(title = "Altitud (m.s.n.m.)")
grid.arrange(h3a, h3b, h3c, ncol = 3)Interpretación: la riqueza de especies (media ≈ 9.7, rango 2-21) y el índice de Shannon (media ≈ 1.76, rango 0.67-2.98) muestran amplia variación entre las 52 parcelas, consistente con el gradiente de intervención antrópica muestreado. La altitud varía entre ~950 y ~1600 m.s.n.m. sin concentraciones extremas.
bb1 <- ggplot(BD_biodiversidad, aes(Habitat, Riqueza, fill = Habitat)) + geom_boxplot(alpha=0.85) +
scale_fill_michi() + theme_michi(9) + theme(axis.text.x = element_text(angle = 25, hjust = 1)) +
labs(title = "Riqueza por hábitat", x = "")
bb2 <- ggplot(BD_biodiversidad, aes(Habitat, Shannon, fill = Habitat)) + geom_boxplot(alpha=0.85) +
scale_fill_michi() + theme_michi(9) + theme(axis.text.x = element_text(angle = 25, hjust = 1)) +
labs(title = "Shannon por hábitat", x = "")
bb3 <- ggplot(BD_biodiversidad, aes(Altitud, Shannon)) + geom_point(color = azul) +
geom_smooth(method = "lm", color = morado, se = FALSE) + theme_michi(9) +
labs(title = "Shannon vs Altitud")
grid.arrange(bb1, bb2, bb3, ncol = 3)cat("Correlación Altitud-Shannon:", round(cor(BD_biodiversidad$Altitud, BD_biodiversidad$Shannon), 3), "\n")## Correlación Altitud-Shannon: 0.752
cat("Correlación Altitud-Riqueza:", round(cor(BD_biodiversidad$Altitud, BD_biodiversidad$Riqueza), 3), "\n")## Correlación Altitud-Riqueza: 0.714
Interpretación: tanto la riqueza como el índice de Shannon decrecen consistentemente a lo largo del gradiente Bosque primario → Bosque secundario → Sistema silvopastoril → Potrero, evidenciando el impacto negativo de la intervención antrópica sobre la biodiversidad de anfibios. Además, existe una correlación positiva moderada-fuerte entre altitud y diversidad de Shannon (r ≈ 0.75): las parcelas más altas del gradiente muestral tienden a coincidir con hábitats menos intervenidos y, por tanto, con mayor diversidad. Es importante señalar que esta correlación debe interpretarse con cautela: la altitud promedio también decrece en el mismo orden que la diversidad de Shannon (Bosque primario ≈1453 m, Bosque secundario ≈1326 m, Sistema silvopastoril ≈1151 m, Potrero ≈1045 m), lo que indica que altitud y tipo de hábitat están confundidos en este diseño observacional — no fueron asignados de forma independiente. Por lo tanto, no puede concluirse causalmente que una mayor altitud por sí misma incremente la diversidad de anfibios; es más plausible que ambas variables reflejen conjuntamente el mismo gradiente de intervención antrópica presente en el muestreo. La diferencia robusta y directamente atribuible es la que se confirma mediante el ANOVA y la prueba LSD entre hábitats, que no depende de este posible factor de confusión.
mod3 <- aov(Shannon ~ Habitat, data = BD_biodiversidad)
tabla_bonita(as.data.frame(summary(mod3)[[1]]), digits = 4, caption = "Tabla ANOVA — Shannon ~ 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 |
res3 <- residuals(mod3)
sw3 <- shapiro.test(res3)
lev3 <- leveneTest(Shannon ~ Habitat, data = BD_biodiversidad)
p_qq3 <- ggplot(data.frame(res3), aes(sample = res3)) + stat_qq(color = azul) + stat_qq_line(color = morado) +
labs(title = "QQ-plot de residuales") + theme_michi(10)
p_res3 <- ggplot(data.frame(fit = fitted(mod3), res = res3), aes(fit, res)) +
geom_point(color = azul) + geom_hline(yintercept = 0, color = morado, linetype = "dashed") +
labs(title = "Residuales vs ajustados") + theme_michi(10)
grid.arrange(p_qq3, p_res3, ncol = 2)Ambos supuestos se satisfacen; el ANOVA y la LSD son válidos sin necesidad de alternativas robustas.
lsd3 <- lsd_test(BD_biodiversidad$Shannon, BD_biodiversidad$Habitat, mod3)
tabla_bonita(lsd3$comparaciones, caption = "Comparaciones LSD por pares — hábitats")| Comparacion | Diferencia | Valor_p | Significativo | |
|---|---|---|---|---|
| Bosque primario | Bosque primario vs Bosque secundario | 0.493 | 0.0001 | Si |
| Bosque primario1 | Bosque primario vs Sistema silvopastoril | 0.865 | 0.0000 | Si |
| Bosque primario2 | Bosque primario vs Potrero | 1.464 | 0.0000 | Si |
| Bosque secundario | Bosque secundario vs Sistema silvopastoril | 0.372 | 0.0018 | Si |
| Bosque secundario1 | Bosque secundario vs Potrero | 0.971 | 0.0000 | Si |
| Sistema silvopastoril | Sistema silvopastoril vs Potrero | 0.599 | 0.0000 | Si |
| Grupo | Media | n | Letra | |
|---|---|---|---|---|
| Bosque primario | Bosque primario | 2.468 | 13 | a |
| Bosque secundario | Bosque secundario | 1.975 | 13 | b |
| Sistema silvopastoril | Sistema silvopastoril | 1.603 | 13 | c |
| Potrero | Potrero | 1.004 | 13 | d |
ggplot(BD_biodiversidad, aes(Habitat, Shannon, fill = Habitat)) +
geom_boxplot(alpha = 0.85) +
scale_fill_michi() +
labs(title = "Diversidad de Shannon-Wiener según tipo de hábitat",
subtitle = "Todos los pares de hábitats difieren significativamente (LSD, p < 0.05)",
x = "", y = "Índice de Shannon-Wiener") +
theme_michi() + theme(axis.text.x = element_text(angle = 15, hjust = 1))Interpretación: el ANOVA detecta diferencias altamente significativas en la diversidad de Shannon entre hábitats (p < 0.001). La prueba LSD muestra que los cuatro hábitats difieren significativamente entre sí, con un patrón decreciente claro: Bosque primario (2.47) > Bosque secundario (1.97) > Sistema silvopastoril (1.60) > Potrero (1.00). Esto confirma que cada nivel de intervención antrópica en el gradiente reduce de forma medible y estadísticamente significativa la diversidad de anfibios, siendo el bosque primario el hábitat que mejor preserva la biodiversidad y el potrero el más degradado en este sentido.