Preparación

Los archivos Salinidad.RData, moluscos.RData y Biodiversidad.RData deben estar en la misma carpeta que este .Rmd.

library(tidyverse)   # manejo de datos y ggplot2
library(car)         # prueba de Levene
library(agricolae)   # prueba LSD de Fisher
library(knitr)       # tablas
library(patchwork)   # combinar gráficos

theme_set(theme_minimal(base_size = 13) +
            theme(plot.title = element_text(face = "bold"),
                  strip.text = element_text(face = "bold")))
colores <- c("#2C7FB8", "#D95F0E", "#31A354", "#756BB1")

load("Salinidad.RData")       # objeto: Salinidad
load("moluscos.RData")        # objeto: BD_moluscos
load("Biodiversidad.RData")   # objeto: BD_biodiversidad

Funciones auxiliares que se usan en los tres puntos (resumen numérico, verificación de supuestos, gráficos de diagnóstico y tamaño de efecto):

# Resumen numérico de una variable
resumen_num <- function(x) {
  tibble(n = sum(!is.na(x)), Media = mean(x), Mediana = median(x),
         DE = sd(x), CV_pct = 100 * sd(x) / mean(x),
         Min = min(x), Max = max(x))
}

# Shapiro-Wilk sobre residuales + Levene
supuestos <- function(modelo, formula, datos) {
  sw <- shapiro.test(residuals(modelo))
  lv <- leveneTest(formula, data = datos)
  tibble(Prueba = c("Shapiro-Wilk (normalidad de residuales)",
                    "Levene (homogeneidad de varianzas)"),
         `Estadístico` = round(c(unname(sw$statistic), lv$`F value`[1]), 3),
         `Valor p` = format.pval(c(sw$p.value, lv$`Pr(>F)`[1]),
                                 digits = 3, eps = 1e-4))
}

# QQ-plot de residuales y residuales vs. ajustados
diagnosticos <- function(modelo) {
  d <- tibble(res = residuals(modelo), ajust = fitted(modelo))
  p1 <- ggplot(d, aes(sample = res)) +
    stat_qq(color = colores[1]) + stat_qq_line(color = colores[2]) +
    labs(title = "QQ-plot de residuales",
         x = "Cuantiles teóricos", y = "Cuantiles muestrales")
  p2 <- ggplot(d, aes(ajust, res)) +
    geom_point(color = colores[1]) +
    geom_hline(yintercept = 0, linetype = "dashed", color = colores[2]) +
    labs(title = "Residuales vs. ajustados",
         x = "Valores ajustados", y = "Residuales")
  p1 + p2
}

# Tamaño de efecto: eta cuadrado por término
eta2 <- function(modelo) {
  tab <- summary(modelo)[[1]]
  ss <- tab$`Sum Sq`
  tibble(`Término` = trimws(rownames(tab)), `Eta²` = round(ss / sum(ss), 3))
}

Punto 1 – Datos Salinidad

Se midieron 45 muestras de suelo en las que crece una planta forrajera natural. La respuesta es la biomasa (g) y las covariables son pH, salinidad, zinc y potasio.

a. Análisis exploratorio univariado

Salinidad %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "valor") %>%
  group_by(Variable) %>%
  summarise(resumen_num(valor), .groups = "drop") %>%
  kable(digits = 2, caption = "Resumen descriptivo de las variables (n = 45)")

Resumen descriptivo de las variables (n = 45)
Variable	n	Media	Mediana	DE	CV_pct	Min	Max
Biomasa	45	1082.17	991.83	546.29	50.48	369.82	2337.33
Potasio	45	797.38	773.30	297.58	37.32	350.73	1441.67
Salinidad	45	30.27	30.00	3.72	12.29	24.00	38.00
Zinc	45	17.83	19.24	8.27	46.40	0.21	31.29
pH	45	4.61	4.45	1.25	27.22	3.20	7.45

Salinidad %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "valor") %>%
  ggplot(aes(valor)) +
  geom_histogram(aes(y = after_stat(density)), bins = 12,
                 fill = colores[1], color = "white", alpha = .85) +
  geom_density(color = colores[2], linewidth = 1) +
  facet_wrap(~Variable, scales = "free", ncol = 2) +
  labs(title = "Distribución de cada característica",
       x = NULL, y = "Densidad")

Salinidad %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "valor") %>%
  ggplot(aes(y = valor)) +
  geom_boxplot(fill = colores[1], alpha = .6, width = .5) +
  facet_wrap(~Variable, scales = "free", nrow = 1) +
  labs(title = "Diagramas de caja", y = NULL) +
  theme(axis.text.x = element_blank())

Interpretación. La biomasa es muy variable (media ≈ 1082 g, DE ≈ 546 g, CV ≈ 50 %) y tiene asimetría positiva: la media supera a la mediana (≈ 992 g), de modo que pocas muestras muy productivas alargan la cola derecha. El pH (media ≈ 4.6) también es asimétrico a la derecha; la mayoría de los suelos son ácidos (mediana 4.45) y solo unos pocos superan pH 6. La salinidad es la más homogénea (CV ≈ 12 %, entre 24 y 38) y aproximadamente simétrica. El zinc (CV ≈ 46 %) muestra asimetría negativa, con unos pocos suelos con concentraciones muy bajas. El potasio varía entre ≈ 350 y 1440 con una ligera asimetría positiva. Que la biomasa sea asimétrica y su dispersión aumente con la media anticipa problemas de homogeneidad de varianzas más adelante.

b. Análisis exploratorio bivariado

Se mide la relación entre la biomasa y pH, salinidad y zinc (el potasio se incluye solo como referencia).

Salinidad %>%
  pivot_longer(c(pH, Salinidad, Zinc, Potasio),
               names_to = "Covariable", values_to = "valor") %>%
  group_by(Covariable) %>%
  summarise(`r de Pearson` = cor(valor, Biomasa),
            `R²` = cor(valor, Biomasa)^2,
            `rho de Spearman` = cor(valor, Biomasa, method = "spearman"),
            .groups = "drop") %>%
  arrange(desc(abs(`r de Pearson`))) %>%
  kable(digits = 3, caption = "Asociación de cada covariable con la biomasa")

Asociación de cada covariable con la biomasa
Covariable	r de Pearson	R²	rho de Spearman
pH	0.928	0.861	0.878
Zinc	-0.781	0.611	-0.625
Potasio	-0.073	0.005	0.102
Salinidad	-0.067	0.004	-0.157

largo <- Salinidad %>%
  select(Biomasa, pH, Salinidad, Zinc) %>%
  pivot_longer(-Biomasa, names_to = "Covariable", values_to = "valor")

cors <- largo %>% group_by(Covariable) %>%
  summarise(r = cor(valor, Biomasa), .groups = "drop")

ggplot(largo, aes(valor, Biomasa)) +
  geom_point(color = colores[1], alpha = .8, size = 2) +
  geom_smooth(method = "lm", color = colores[2], fill = "#FDAE6B") +
  geom_text(data = cors, inherit.aes = FALSE,
            aes(x = -Inf, y = Inf, label = paste0("r = ", round(r, 2))),
            hjust = -.2, vjust = 1.6, fontface = "bold") +
  facet_wrap(~Covariable, scales = "free_x") +
  labs(title = "Biomasa vs. covariables del suelo",
       x = NULL, y = "Biomasa (g)")

cor(Salinidad$pH, Salinidad$Zinc)   # relación entre las dos covariables más asociadas

## [1] -0.7204699

Interpretación. El pH es la covariable que presenta mayor relación con la biomasa (r ≈ +0.93; explica ≈ 86 % de su variabilidad lineal): suelos menos ácidos se asocian con más biomasa. Le sigue el zinc (r ≈ −0.78), con relación negativa: más zinc, menos biomasa. La salinidad prácticamente no se relaciona con la biomasa (r ≈ −0.07) en el rango observado. Como pH y zinc están correlacionados entre sí (r ≈ −0.72), parte del efecto aparente del zinc puede deberse al pH (en suelos más ácidos el zinc suele ser más soluble y disponible), por lo que no deben interpretarse como efectos independientes. Por lo tanto, la variable a categorizar en el literal (c) es el pH.

c. ANOVA de una vía según el nivel de pH

Categorización en terciles

Se divide el pH en tres niveles con los terciles (cada nivel con la misma cantidad de muestras, 15 por grupo), de forma que la comparación tenga tamaños iguales.

cortes <- quantile(Salinidad$pH, probs = c(0, 1/3, 2/3, 1))
Salinidad <- Salinidad %>%
  mutate(pH_cat = cut(pH, breaks = cortes, include.lowest = TRUE,
                      labels = c("Bajo", "Medio", "Alto")),
         logBiomasa = log(Biomasa))

Salinidad %>%
  group_by(pH_cat) %>%
  summarise(n = n(), `pH mín` = min(pH), `pH máx` = max(pH),
            Media = mean(Biomasa), Mediana = median(Biomasa),
            DE = sd(Biomasa), CV_pct = 100 * sd(Biomasa) / mean(Biomasa)) %>%
  kable(digits = 2, caption = "Biomasa (g) por nivel de pH")

Biomasa (g) por nivel de pH
pH_cat	n	pH mín	pH máx	Media	Mediana	DE	CV_pct
Bajo	15	3.20	3.75	593.02	545.54	164.83	27.79
Medio	15	3.95	4.85	1048.11	1039.64	244.91	23.37
Alto	15	5.00	7.45	1605.39	1422.84	547.60	34.11

ggplot(Salinidad, aes(pH_cat, Biomasa, fill = pH_cat)) +
  geom_boxplot(alpha = .7, outlier.shape = NA) +
  geom_jitter(width = .1, alpha = .6) +
  scale_fill_manual(values = colores[c(2, 3, 1)], guide = "none") +
  labs(title = "Biomasa según nivel de pH", x = "Nivel de pH", y = "Biomasa (g)")

Verificación de supuestos del modelo original

modelo_raw <- aov(Biomasa ~ pH_cat, data = Salinidad)
summary(modelo_raw)

##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## pH_cat       2 7712683 3856342   29.89 8.45e-09 ***
## Residuals   42 5418235  129006                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

supuestos(modelo_raw, Biomasa ~ pH_cat, Salinidad) %>% kable()

Prueba	Estadístico	Valor p
Shapiro-Wilk (normalidad de residuales)	0.973	0.374885
Levene (homogeneidad de varianzas)	8.675	0.000702

diagnosticos(modelo_raw)

Los residuales son compatibles con la normalidad, pero la prueba de Levene rechaza la homogeneidad de varianzas (la desviación estándar pasa de ≈ 165 g en pH bajo a ≈ 548 g en pH alto). Como el ANOVA clásico asume varianzas iguales, no se interpreta el modelo sobre la biomasa original. Se aplica una transformación logarítmica, que estabiliza la varianza cuando la dispersión crece con la media.

Modelo con la biomasa transformada (log)

modelo_log <- aov(logBiomasa ~ pH_cat, data = Salinidad)
summary(modelo_log)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## pH_cat       2  7.154   3.577   40.81 1.42e-10 ***
## Residuals   42  3.681   0.088                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

eta2(modelo_log) %>% kable(caption = "Tamaño de efecto")

Tamaño de efecto
Término	Eta²
pH_cat	0.66
Residuals	0.34

supuestos(modelo_log, logBiomasa ~ pH_cat, Salinidad) %>% kable()

Prueba	Estadístico	Valor p
Shapiro-Wilk (normalidad de residuales)	0.979	0.578
Levene (homogeneidad de varianzas)	1.413	0.255

diagnosticos(modelo_log)

Con log(Biomasa) se cumplen ambos supuestos (Shapiro-Wilk y Levene no rechazan), por lo que el ANOVA es válido. Como verificación adicional, una prueba no paramétrica que no necesita estos supuestos:

kruskal.test(Biomasa ~ pH_cat, data = Salinidad)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Biomasa by pH_cat
## Kruskal-Wallis chi-squared = 29.2, df = 2, p-value = 4.564e-07

Comparaciones post-hoc (LSD de Fisher)

lsd1 <- LSD.test(modelo_log, "pH_cat", console = TRUE)

## 
## Study: modelo_log ~ "pH_cat"
## 
## LSD t Test for logBiomasa 
## 
## Mean Square Error:  0.08763569 
## 
## pH_cat,  means and individual ( 95 %) CI
## 
##       logBiomasa       std  r         se      LCL      UCL      Min      Max
## Alto    7.322857 0.3608580 15 0.07643546 7.168604 7.477110 6.640242 7.756763
## Bajo    6.351728 0.2641536 15 0.07643546 6.197475 6.505981 5.913025 6.885014
## Medio   6.926945 0.2508216 15 0.07643546 6.772692 7.081198 6.342922 7.307387
##            Q25      Q50      Q75
## Alto  7.089255 7.260407 7.692731
## Bajo  6.176084 6.301772 6.491777
## Medio 6.789658 6.946627 7.088730
## 
## Alpha: 0.05 ; DF Error: 42
## Critical Value of t: 2.018082 
## 
## least Significant Difference: 0.2181467 
## 
## Treatments with the same letter are not significantly different.
## 
##       logBiomasa groups
## Alto    7.322857      a
## Medio   6.926945      b
## Bajo    6.351728      c

letras1 <- as.data.frame(lsd1$groups) %>%
  rownames_to_column("pH_cat") %>%
  left_join(Salinidad %>% group_by(pH_cat) %>%
              summarise(ymax = max(Biomasa), .groups = "drop"),
            by = "pH_cat")

Salinidad %>%
  group_by(pH_cat) %>%
  summarise(`Media geométrica (g)` = exp(mean(logBiomasa))) %>%
  kable(digits = 1, caption = "Media geométrica de la biomasa (equivale a la media de log(Biomasa) en la escala original)")

Media geométrica de la biomasa (equivale a la media de log(Biomasa) en la escala original)
pH_cat	Media geométrica (g)
Bajo	573.5
Medio	1019.4
Alto	1514.5

ggplot(Salinidad, aes(pH_cat, Biomasa, fill = pH_cat)) +
  geom_boxplot(alpha = .7, outlier.shape = NA) +
  geom_jitter(width = .1, alpha = .5) +
  geom_text(data = letras1, aes(y = ymax * 1.06, label = groups),
            fontface = "bold", size = 6) +
  scale_fill_manual(values = colores[c(2, 3, 1)], guide = "none") +
  labs(title = "Biomasa por nivel de pH (letras distintas = diferencia significativa, LSD)",
       x = "Nivel de pH", y = "Biomasa (g)")

Conclusión del punto 1. El nivel de pH del suelo tiene un efecto altamente significativo sobre la biomasa (ANOVA sobre log-biomasa, p < 0.001; el nivel de pH explica ≈ 66 % de la variación, un efecto muy grande). El LSD muestra que los tres niveles difieren entre sí: la biomasa aumenta de forma escalonada de pH bajo (≈ 593 g) a medio (≈ 1048 g) y a alto (≈ 1605 g). En términos del problema, la planta forrajera produce mucho más en suelos menos ácidos, y pasar de un suelo muy ácido a uno moderadamente ácido ya casi duplica la biomasa. Esto concuerda con el análisis bivariado y es relevante para el manejo del suelo (por ejemplo, encalado en suelos muy ácidos). Al ser un estudio observacional, la relación es una asociación y no prueba causalidad.

Punto 2 – Datos Moluscos

Dos tipos de moluscos (A y B) se sometieron a tres concentraciones de agua de mar (100 %, 75 % y 50 %) y se midió el consumo de oxígeno por unidad de peso seco.

moluscos <- BD_moluscos %>%
  mutate(c_agua  = factor(c_agua, levels = c(100, 75, 50),
                          labels = c("100%", "75%", "50%")),
         molusco = factor(molusco))

a. Análisis exploratorio univariado

moluscos %>% summarise(resumen_num(cons_o)) %>%
  kable(digits = 2, caption = "Consumo de oxígeno (todas las observaciones)")

Consumo de oxígeno (todas las observaciones)
n	Media	Mediana	DE	CV_pct	Min	Max
48	9.3	9.7	3.68	39.58	1.8	18.8

moluscos %>% count(molusco, c_agua) %>%
  pivot_wider(names_from = c_agua, values_from = n) %>%
  kable(caption = "Número de observaciones por combinación (diseño balanceado)")

Número de observaciones por combinación (diseño balanceado)
molusco	100%	75%	50%
A	8	8	8
B	8	8	8

p1 <- ggplot(moluscos, aes(cons_o)) +
  geom_histogram(aes(y = after_stat(density)), bins = 10,
                 fill = colores[1], color = "white", alpha = .85) +
  geom_density(color = colores[2], linewidth = 1) +
  labs(title = "Consumo de oxígeno", x = "Consumo de O2 / peso seco", y = "Densidad")
p2 <- ggplot(moluscos, aes(y = cons_o)) +
  geom_boxplot(fill = colores[1], alpha = .6, width = .4) +
  labs(title = "Diagrama de caja", y = NULL) +
  theme(axis.text.x = element_blank())
p1 + p2 + plot_layout(widths = c(2, 1))

Interpretación. El consumo de oxígeno tiene una media de ≈ 9.3 y mediana ≈ 9.7 (DE ≈ 3.7, CV ≈ 40 %), con una distribución aproximadamente simétrica entre 1.8 y 18.8. Hay variabilidad considerable entre individuos, y los valores altos (> 14) son poco frecuentes pero no extremos. El diseño es balanceado: 8 individuos por cada combinación de molusco y concentración (48 en total), lo que facilita el análisis.

b. Análisis exploratorio bivariado

moluscos %>%
  group_by(molusco, c_agua) %>%
  summarise(n = n(), Media = mean(cons_o), Mediana = median(cons_o),
            DE = sd(cons_o), .groups = "drop") %>%
  kable(digits = 2, caption = "Consumo de oxígeno por tipo de molusco y concentración")

Consumo de oxígeno por tipo de molusco y concentración
molusco	c_agua	n	Media	Mediana	DE
A	100%	8	9.94	9.30	2.75
A	75%	8	7.89	7.18	2.74
A	50%	8	12.18	11.11	3.09
B	100%	8	7.41	6.14	2.84
B	75%	8	6.10	5.60	2.74
B	50%	8	12.33	12.85	3.52

moluscos %>% group_by(c_agua) %>%
  summarise(Media = mean(cons_o), DE = sd(cons_o)) %>%
  kable(digits = 2, caption = "Por concentración (ambos moluscos)")

Por concentración (ambos moluscos)
c_agua	Media	DE
100%	8.67	3.0
75%	6.99	2.8
50%	12.25	3.2

g1 <- ggplot(moluscos, aes(c_agua, cons_o, fill = molusco)) +
  geom_boxplot(alpha = .75, position = position_dodge(.8)) +
  scale_fill_manual(values = colores[1:2], name = "Molusco") +
  labs(title = "Consumo de O2 por concentración y tipo",
       x = "Concentración de agua de mar", y = "Consumo de O2 / peso seco")

g2 <- moluscos %>%
  group_by(molusco, c_agua) %>%
  summarise(m = mean(cons_o), se = sd(cons_o) / sqrt(n()), .groups = "drop") %>%
  ggplot(aes(c_agua, m, color = molusco, group = molusco)) +
  geom_line(linewidth = 1) + geom_point(size = 3) +
  geom_errorbar(aes(ymin = m - se, ymax = m + se), width = .12) +
  scale_color_manual(values = colores[1:2], name = "Molusco") +
  labs(title = "Gráfico de interacción (media ± EE)",
       x = "Concentración de agua de mar", y = "Consumo medio de O2")
g1 + g2

Interpretación. El consumo de oxígeno no cambia de forma lineal con la concentración: es mayor al 50 % (≈ 12.3), menor al 75 % (≈ 7.0) y intermedio al 100 % (≈ 8.7). Las conclusiones son similares para ambos moluscos: las dos curvas siguen el mismo patrón en forma de “V” (alto en 50 %, mínimo en 75 %), con líneas casi paralelas. El molusco A consume algo más que el B en 100 % y 75 %, pero en 50 % ambos consumen prácticamente lo mismo (≈ 12.2 vs. 12.3). La dispersión dentro de cada grupo es considerable, de modo que hace falta el ANOVA para juzgar qué diferencias son reales.

c. ANOVA de dos vías con interacción

modelo2 <- aov(cons_o ~ molusco * c_agua, data = moluscos)

Verificación de supuestos

supuestos(modelo2, cons_o ~ molusco * c_agua, moluscos) %>% kable()

Prueba	Estadístico	Valor p
Shapiro-Wilk (normalidad de residuales)	0.958	0.0857
Levene (homogeneidad de varianzas)	0.172	0.9715

diagnosticos(modelo2)

Los residuales no rechazan normalidad (el valor p de Shapiro-Wilk queda por encima de 0.05) y Levene no rechaza la homogeneidad de varianzas (p ≈ 0.97). Se cumplen los supuestos y el ANOVA clásico es apropiado.

Tabla ANOVA

summary(modelo2)

##                Df Sum Sq Mean Sq F value   Pr(>F)    
## molusco         1   23.2   23.23   2.651    0.111    
## c_agua          2  230.8  115.41  13.171 3.63e-05 ***
## molusco:c_agua  2   15.4    7.68   0.876    0.424    
## Residuals      42  368.0    8.76                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

eta2(modelo2) %>% kable(caption = "Tamaño de efecto (eta²)")

Tamaño de efecto (eta²)
Término	Eta²
molusco	0.036
c_agua	0.362
molusco:c_agua	0.024
Residuals	0.577

Comparaciones post-hoc (LSD)

Como la interacción no es significativa, la lectura principal se hace sobre los efectos principales; se hace LSD para la concentración de agua (el único factor significativo). Como complemento se presenta el LSD entre las seis combinaciones.

lsd2 <- LSD.test(modelo2, "c_agua", console = TRUE)

## 
## Study: modelo2 ~ "c_agua"
## 
## LSD t Test for cons_o 
## 
## Mean Square Error:  8.762171 
## 
## c_agua,  means and individual ( 95 %) CI
## 
##        cons_o      std  r        se       LCL       UCL  Min  Max    Q25    Q50
## 100%  8.67125 3.000940 16 0.7400241  7.177821 10.164679 3.68 14.0  6.140  8.595
## 50%  12.25062 3.199643 16 0.7400241 10.757196 13.744054 6.38 18.8 10.085 11.455
## 75%   6.99250 2.804093 16 0.7400241  5.499071  8.485929 1.80 13.2  5.200  6.430
##          Q75
## 100% 10.5750
## 50%  14.5000
## 75%   8.7675
## 
## Alpha: 0.05 ; DF Error: 42
## Critical Value of t: 2.018082 
## 
## least Significant Difference: 2.112028 
## 
## Treatments with the same letter are not significantly different.
## 
##        cons_o groups
## 50%  12.25062      a
## 100%  8.67125      b
## 75%   6.99250      b

moluscos <- moluscos %>% mutate(tratamiento = interaction(molusco, c_agua))
modelo2_trat <- aov(cons_o ~ tratamiento, data = moluscos)
lsd2_trat <- LSD.test(modelo2_trat, "tratamiento", console = TRUE)

## 
## Study: modelo2_trat ~ "tratamiento"
## 
## LSD t Test for cons_o 
## 
## Mean Square Error:  8.762171 
## 
## tratamiento,  means and individual ( 95 %) CI
## 
##          cons_o      std r       se       LCL       UCL  Min   Max     Q25
## A.100%  9.93625 2.747976 8 1.046552  7.824222 12.048278 6.78 14.00  7.9850
## A.50%  12.17500 3.090178 8 1.046552 10.062972 14.287028 9.74 18.80 10.3100
## A.75%   7.89000 2.739578 8 1.046552  5.777972 10.002028 5.20 13.20  6.0775
## B.100%  7.40625 2.844076 8 1.046552  5.294222  9.518278 3.68 11.60  5.7225
## B.50%  12.32625 3.517909 8 1.046552 10.214222 14.438278 6.38 17.70 10.0575
## B.75%   6.09500 2.739108 8 1.046552  3.982972  8.207028 1.80  9.96  4.8300
##           Q50     Q75
## A.100%  9.295 11.7250
## A.50%  11.110 12.5000
## A.75%   7.180  8.8925
## B.100%  6.140 10.1000
## B.50%  12.850 14.5000
## B.75%   5.595  7.3425
## 
## Alpha: 0.05 ; DF Error: 42
## Critical Value of t: 2.018082 
## 
## least Significant Difference: 2.986858 
## 
## Treatments with the same letter are not significantly different.
## 
##          cons_o groups
## B.50%  12.32625      a
## A.50%  12.17500      a
## A.100%  9.93625     ab
## A.75%   7.89000     bc
## B.100%  7.40625     bc
## B.75%   6.09500      c

letras2 <- as.data.frame(lsd2$groups) %>%
  rownames_to_column("c_agua") %>%
  left_join(moluscos %>% group_by(c_agua) %>%
              summarise(ymax = max(cons_o), .groups = "drop"), by = "c_agua")

ggplot(moluscos, aes(c_agua, cons_o, fill = c_agua)) +
  geom_boxplot(alpha = .7, outlier.shape = NA) +
  geom_jitter(aes(shape = molusco), width = .12, alpha = .7, size = 2) +
  geom_text(data = letras2, aes(y = ymax * 1.07, label = groups),
            fontface = "bold", size = 6) +
  scale_fill_manual(values = colores[c(1, 3, 2)], guide = "none") +
  labs(title = "Consumo de O2 por concentración (letras distintas = diferencia, LSD)",
       x = "Concentración de agua de mar", y = "Consumo de O2 / peso seco",
       shape = "Molusco")

Conclusión del punto 2. Interacción (molusco × concentración): no es significativa (p ≈ 0.42), es decir, el efecto de la concentración de agua de mar es el mismo para los dos tipos de molusco y viceversa. Tipo de molusco: no tiene un efecto principal significativo (p ≈ 0.11); A consume en promedio algo más que B (≈ 10.0 vs. 8.6) pero la diferencia no se puede distinguir de la variación natural. Concentración de agua de mar: efecto altamente significativo (p < 0.001), explica ≈ 36 % de la variación total. Según el LSD, el consumo al 50 % es mayor que al 100 % y al 75 %, mientras que 100 % y 75 % no difieren entre sí. Biológicamente, el mayor consumo de oxígeno en agua más diluida es consistente con un mayor costo energético para mantener el equilibrio osmótico (osmorregulación) cuando el medio es hipotónico, y esto ocurre de forma similar en ambos moluscos. El mínimo observado en 75 % (menor que en 100 %) no es estadísticamente distinguible del de 100 %, así que no debe interpretarse como una tendencia real.

Punto 3 – Datos Biodiversidad

Se evalúa el efecto del uso del suelo sobre la biodiversidad de anfibios en una reserva forestal: 52 parcelas (13 por hábitat) en un gradiente de intervención antrópica.

bio <- BD_biodiversidad %>%
  mutate(Habitat = factor(Habitat,
                          levels = c("Bosque primario", "Bosque secundario",
                                     "Sistema silvopastoril", "Potrero")))

a. Análisis exploratorio univariado

bio %>%
  select(Riqueza, Shannon, Altitud) %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "valor") %>%
  group_by(Variable) %>%
  summarise(resumen_num(valor), .groups = "drop") %>%
  kable(digits = 2, caption = "Resumen descriptivo (52 parcelas)")

Resumen descriptivo (52 parcelas)
Variable	n	Media	Mediana	DE	CV_pct	Min	Max
Altitud	52	1243.89	1254.80	172.71	13.88	950.90	1600.80
Riqueza	52	9.73	9.50	4.37	44.92	2.00	21.00
Shannon	52	1.76	1.85	0.61	34.44	0.67	2.98

bio %>%
  select(Riqueza, Shannon, Altitud) %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "valor") %>%
  ggplot(aes(valor)) +
  geom_histogram(aes(y = after_stat(density)), bins = 12,
                 fill = colores[3], color = "white", alpha = .85) +
  geom_density(color = colores[2], linewidth = 1) +
  facet_wrap(~Variable, scales = "free", ncol = 2) +
  labs(title = "Distribución de riqueza, diversidad de Shannon y altitud",
       x = NULL, y = "Densidad")

Interpretación. La riqueza promedio es ≈ 9.7 especies por parcela (DE ≈ 4.4), con un máximo de 21. El índice de Shannon promedia ≈ 1.76 (DE ≈ 0.61; rango 0.67 a 2.98), y la altitud ≈ 1244 m (DE ≈ 173 m; rango ≈ 951 a 1601 m). Estos valores mezclan cuatro hábitats con condiciones muy distintas, y ese contraste se explora en el literal (b).

b. Análisis exploratorio bivariado

bio %>%
  group_by(Habitat) %>%
  summarise(n = n(),
            `Riqueza (media)` = mean(Riqueza), `Riqueza (DE)` = sd(Riqueza),
            `Shannon (media)` = mean(Shannon), `Shannon (DE)` = sd(Shannon),
            `Altitud (media)` = mean(Altitud), `Altitud (DE)` = sd(Altitud)) %>%
  kable(digits = 2, caption = "Resumen por tipo de hábitat")

Resumen por tipo de hábitat
Habitat	n	Riqueza (media)	Riqueza (DE)	Shannon (media)	Shannon (DE)	Altitud (media)	Altitud (DE)
Bosque primario	13	14.77	2.74	2.47	0.38	1453.19	77.50
Bosque secundario	13	11.08	3.01	1.97	0.22	1326.15	40.60
Sistema silvopastoril	13	8.08	2.50	1.60	0.28	1151.48	97.89
Potrero	13	5.00	1.29	1.00	0.24	1044.75	50.32

b1 <- ggplot(bio, aes(Habitat, Riqueza, fill = Habitat)) +
  geom_boxplot(alpha = .75) + geom_jitter(width = .1, alpha = .5) +
  scale_fill_manual(values = colores, guide = "none") +
  labs(title = "Riqueza de especies", x = NULL, y = "Número de especies") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
b2 <- ggplot(bio, aes(Habitat, Shannon, fill = Habitat)) +
  geom_boxplot(alpha = .75) + geom_jitter(width = .1, alpha = .5) +
  scale_fill_manual(values = colores, guide = "none") +
  labs(title = "Diversidad de Shannon", x = NULL, y = "Índice de Shannon") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
b1 + b2

ggplot(bio, aes(Altitud, Shannon)) +
  geom_point(aes(color = Habitat), size = 2.5, alpha = .85) +
  geom_smooth(method = "lm", color = "black", se = FALSE, linewidth = 1) +
  geom_smooth(aes(color = Habitat), method = "lm", se = FALSE,
              linewidth = .6, linetype = "dashed") +
  scale_color_manual(values = colores) +
  labs(title = "Altitud vs. diversidad de Shannon",
       subtitle = "Línea negra: todas las parcelas; líneas punteadas: dentro de cada hábitat",
       x = "Altitud (m.s.n.m.)", y = "Índice de Shannon")

cat("Correlación global Altitud-Shannon:",
    round(cor(bio$Altitud, bio$Shannon), 3), "\n")

## Correlación global Altitud-Shannon: 0.752

bio %>% group_by(Habitat) %>%
  summarise(`r Altitud-Shannon` = cor(Altitud, Shannon), n = n()) %>%
  kable(digits = 2, caption = "Correlación altitud-Shannon dentro de cada hábitat")

Correlación altitud-Shannon dentro de cada hábitat
Habitat	r Altitud-Shannon	n
Bosque primario	-0.33	13
Bosque secundario	-0.01	13
Sistema silvopastoril	-0.51	13
Potrero	0.03	13

Interpretación. Tanto la riqueza como la diversidad de Shannon disminuyen de forma ordenada a lo largo del gradiente de intervención: bosque primario (riqueza ≈ 14.8; Shannon ≈ 2.47) > bosque secundario (≈ 11.1; 1.97) > sistema silvopastoril (≈ 8.1; 1.60) > potrero (≈ 5.0; 1.00). La dispersión dentro de cada hábitat es moderada y los hábitats se separan bastante bien. Respecto a la altitud, la correlación global con Shannon es alta (r ≈ 0.75), pero hay que ser cuidadoso: los hábitats también difieren en altitud (bosque primario ≈ 1453 m, potrero ≈ 1045 m), de modo que altitud y hábitat están confundidos. Al mirar dentro de cada hábitat la relación desaparece o incluso se vuelve negativa (r entre ≈ −0.5 y 0.03, con solo 13 parcelas por hábitat). Es decir, la asociación global se debe sobre todo a que los hábitats más conservados están ubicados a mayor altitud, y no a un efecto de la altitud en sí.

c. ANOVA de una vía: Shannon según hábitat

modelo3 <- aov(Shannon ~ Habitat, data = bio)

Verificación de supuestos

supuestos(modelo3, Shannon ~ Habitat, bio) %>% kable()

Prueba	Estadístico	Valor p
Shapiro-Wilk (normalidad de residuales)	0.978	0.463
Levene (homogeneidad de varianzas)	1.277	0.293

diagnosticos(modelo3)

Ni Shapiro-Wilk ni Levene rechazan sus hipótesis nulas (p > 0.05): los residuales son compatibles con normalidad y las varianzas son homogéneas entre hábitats. El ANOVA clásico es apropiado.

Tabla ANOVA y tamaño de efecto

summary(modelo3)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Habitat      3 14.862   4.954   60.61 2.38e-16 ***
## Residuals   48  3.923   0.082                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

eta2(modelo3) %>% kable(caption = "Tamaño de efecto (eta²)")

Tamaño de efecto (eta²)
Término	Eta²
Habitat	0.791
Residuals	0.209

Comparaciones post-hoc (LSD de Fisher)

El ANOVA es significativo, por lo que se justifican las comparaciones múltiples.

lsd3 <- LSD.test(modelo3, "Habitat", console = TRUE)

## 
## Study: modelo3 ~ "Habitat"
## 
## LSD t Test for Shannon 
## 
## Mean Square Error:  0.08173622 
## 
## Habitat,  means and individual ( 95 %) CI
## 
##                        Shannon       std  r         se       LCL      UCL  Min
## Bosque primario       2.467692 0.3765668 13 0.07929314 2.3082628 2.627122 1.73
## Bosque secundario     1.974615 0.2204395 13 0.07929314 1.8151858 2.134045 1.44
## Potrero               1.003846 0.2396712 13 0.07929314 0.8444166 1.163276 0.67
## Sistema silvopastoril 1.603077 0.2812586 13 0.07929314 1.4436474 1.762506 1.14
##                        Max  Q25  Q50  Q75
## Bosque primario       2.98 2.13 2.53 2.75
## Bosque secundario     2.27 1.86 2.03 2.10
## Potrero               1.34 0.80 0.98 1.20
## Sistema silvopastoril 2.06 1.41 1.62 1.86
## 
## Alpha: 0.05 ; DF Error: 48
## Critical Value of t: 2.010635 
## 
## least Significant Difference: 0.2254674 
## 
## Treatments with the same letter are not significantly different.
## 
##                        Shannon groups
## Bosque primario       2.467692      a
## Bosque secundario     1.974615      b
## Sistema silvopastoril 1.603077      c
## Potrero               1.003846      d

letras3 <- as.data.frame(lsd3$groups) %>%
  rownames_to_column("Habitat") %>%
  left_join(bio %>% group_by(Habitat) %>%
              summarise(ymax = max(Shannon), .groups = "drop"), by = "Habitat") %>%
  mutate(Habitat = factor(Habitat, levels = levels(bio$Habitat)))

ggplot(bio, aes(Habitat, Shannon, fill = Habitat)) +
  geom_boxplot(alpha = .75, outlier.shape = NA) +
  geom_jitter(width = .1, alpha = .5) +
  geom_text(data = letras3, aes(y = ymax + .15, label = groups),
            fontface = "bold", size = 6) +
  scale_fill_manual(values = colores, guide = "none") +
  labs(title = "Diversidad de Shannon por hábitat (letras distintas = diferencia, LSD)",
       x = NULL, y = "Índice de Shannon") +
  theme(axis.text.x = element_text(angle = 15, hjust = 1))

Conclusión del punto 3. El tipo de hábitat tiene un efecto muy fuerte y significativo sobre la diversidad de Shannon (F ≈ 60.6, p < 0.001; el hábitat explica ≈ 79 % de la variación). El LSD indica que los cuatro hábitats difieren entre sí: bosque primario (≈ 2.47) > bosque secundario (≈ 1.97) > sistema silvopastoril (≈ 1.60) > potrero (≈ 1.00). La diversidad de anfibios cae de forma escalonada a medida que aumenta la intervención. El potrero conserva menos de la mitad de la diversidad del bosque primario, y el bosque secundario aún no recupera el nivel del primario. El sistema silvopastoril se comporta como un punto intermedio: mantiene más diversidad que el potrero (diferencia ≈ 0.6), lo que lo hace una opción de uso productivo menos dañina para los anfibios, aunque sin llegar al nivel del bosque. Por tratarse de un estudio observacional, y dado que la altitud difiere entre hábitats, el efecto del hábitat no puede separarse por completo del de la altitud.

Conclusiones generales

Punto 1: el pH del suelo es el mejor predictor de la biomasa de la planta forrajera y su nivel genera tres grupos distintos de producción. Fue necesario transformar la respuesta (log) para cumplir el supuesto de varianzas homogéneas.
Punto 2: la concentración de agua de mar afecta el consumo de oxígeno, igual en ambos moluscos (sin interacción); el mayor consumo ocurre en la menor concentración (50 %).
Punto 3: la diversidad de anfibios disminuye de forma escalonada con la intervención antrópica, y los cuatro hábitats difieren entre sí.

Información de la sesión

sessionInfo()

## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.utf8    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] patchwork_1.3.2 knitr_1.52      agricolae_1.3-7 car_3.1-5      
##  [5] carData_3.0-6   lubridate_1.9.5 forcats_1.0.1   stringr_1.6.0  
##  [9] dplyr_1.2.1     purrr_1.2.2     readr_2.1.6     tidyr_1.3.2    
## [13] tibble_3.3.1    ggplot2_4.0.3   tidyverse_2.0.0
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        generics_0.1.4     stringi_1.8.7      lattice_0.22-7    
##  [5] hms_1.1.4          digest_0.6.39      magrittr_2.0.4     evaluate_1.0.5    
##  [9] grid_4.5.2         timechange_0.4.0   RColorBrewer_1.1-3 fastmap_1.2.0     
## [13] Matrix_1.7-4       jsonlite_2.0.0     Formula_1.2-6      mgcv_1.9-3        
## [17] scales_1.4.0       jquerylib_0.1.4    abind_1.4-8        cli_3.6.5         
## [21] rlang_1.1.7        splines_4.5.2      withr_3.0.2        cachem_1.1.0      
## [25] yaml_2.3.12        otel_0.2.0         tools_4.5.2        AlgDesign_1.2.1.2 
## [29] tzdb_0.5.0         vctrs_0.7.1        R6_2.6.1           lifecycle_1.0.5   
## [33] MASS_7.3-65        cluster_2.1.8.1    pkgconfig_2.0.3    pillar_1.11.1     
## [37] bslib_0.10.0       gtable_0.3.6       glue_1.8.0         xfun_0.56         
## [41] tidyselect_1.2.1   rstudioapi_0.19.0  farver_2.1.2       htmltools_0.5.9   
## [45] nlme_3.1-168       labeling_0.4.3     rmarkdown_2.30     compiler_4.5.2    
## [49] S7_0.2.2

Taller 1: Análisis exploratorio, ANOVA y pruebas post-ANOVA

Sebastian Alexis Lopez Cespedes y Santiago Ramirez Lopez

03/10/2026

Preparación

Punto 1 – Datos Salinidad

a. Análisis exploratorio univariado

b. Análisis exploratorio bivariado

c. ANOVA de una vía según el nivel de pH

Categorización en terciles

Verificación de supuestos del modelo original

Modelo con la biomasa transformada (log)

Comparaciones post-hoc (LSD de Fisher)

Punto 2 – Datos Moluscos

a. Análisis exploratorio univariado

b. Análisis exploratorio bivariado

c. ANOVA de dos vías con interacción

Verificación de supuestos

Tabla ANOVA

Comparaciones post-hoc (LSD)

Punto 3 – Datos Biodiversidad

a. Análisis exploratorio univariado

b. Análisis exploratorio bivariado

c. ANOVA de una vía: Shannon según hábitat

Verificación de supuestos

Tabla ANOVA y tamaño de efecto

Comparaciones post-hoc (LSD de Fisher)

Conclusiones generales

Información de la sesión