1 Paquetes y carga de datos

# install.packages(c("agricolae","car","ggplot2","dplyr","tidyr","moments","knitr"))
library(agricolae)   # prueba LSD
library(car)         # prueba de Levene
library(ggplot2)
library(dplyr)
library(tidyr)
library(moments)     # asimetría y curtosis
library(knitr)

La base de datos contiene 45 observaciones y 5 variables numéricas, sin datos faltantes (0 valores ausentes en total): Biomasa (respuesta, en gramos) y las covariables del suelo pH, Sal_suelo, Zinc y Potasio.

Primeras filas de la base de datos
Biomasa pH Sal_suelo 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 Punto (a) — Análisis exploratorio univariado

El objetivo es: describir cada variable por separado, sin relacionarla con ninguna otra. Interesan cuatro cosas: dónde se ubica el centro (media, mediana), qué tan dispersos están los datos (desviación estándar, coeficiente de variación, rango intercuartílico), qué forma tiene la distribución (simetría) y si existen valores atípicos.

2.1 Tabla de estadísticas descriptivas

resumen <- data.frame(
  Variable = names(Salinidad),
  Media    = sapply(Salinidad, mean),
  Mediana  = sapply(Salinidad, median),
  DE       = sapply(Salinidad, sd),
  CV_pct   = sapply(Salinidad, function(x) 100 * sd(x) / mean(x)),
  Minimo   = sapply(Salinidad, min),
  Q1       = sapply(Salinidad, function(x) quantile(x, 0.25)),
  Q3       = sapply(Salinidad, function(x) quantile(x, 0.75)),
  Maximo   = sapply(Salinidad, max),
  Asimetria= sapply(Salinidad, skewness),
  row.names = NULL
)
kable(resumen, digits = 3,
      caption = "Estadísticas descriptivas de las 45 muestras")
Estadísticas descriptivas de las 45 muestras
Variable Media Mediana DE CV_pct Minimo Q1 Q3 Maximo Asimetria
Biomasa 1082.173 991.829 546.287 50.481 369.823 654.825 1346.880 2337.326 0.937
pH 4.609 4.450 1.255 27.224 3.200 3.450 5.350 7.450 0.904
Sal_suelo 30.267 30.000 3.720 12.290 24.000 27.000 33.000 38.000 0.318
Zinc 17.831 19.242 8.274 46.404 0.210 13.985 22.676 31.286 -0.686
Potasio 797.378 773.300 297.576 37.319 350.730 526.970 954.110 1441.670 0.496

2.2 Histogramas y densidades

Salinidad_largo <- pivot_longer(Salinidad, everything(),
                                names_to = "Variable", values_to = "Valor")

ggplot(Salinidad_largo, aes(x = Valor)) +
  geom_histogram(aes(y = after_stat(density)), bins = 10,
                 fill = "steelblue", color = "white") +
  geom_density(color = "darkred", linewidth = 0.9) +
  facet_wrap(~ Variable, scales = "free", ncol = 3) +
  labs(title = "Distribución univariada de cada variable",
       x = NULL, y = "Densidad") +
  theme_bw()

2.3 Diagramas de caja

ggplot(Salinidad_largo, aes(y = Valor)) +
  geom_boxplot(fill = "lightsteelblue", outlier.color = "red") +
  facet_wrap(~ Variable, scales = "free_y", ncol = 5) +
  labs(title = "Diagramas de caja por variable", y = NULL) +
  theme_bw() +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

2.4 Prueba de normalidad por variable (complementaria)

shapiro_tabla <- data.frame(
  Variable = names(Salinidad),
  p_valor  = sapply(Salinidad, function(x) shapiro.test(x)$p.value)
)
kable(shapiro_tabla, digits = 4,
      caption = "Shapiro-Wilk por variable (H0: la variable es normal)")
Shapiro-Wilk por variable (H0: la variable es normal)
Variable p_valor
Biomasa Biomasa 0.0008
pH pH 0.0003
Sal_suelo Sal_suelo 0.0507
Zinc Zinc 0.0091
Potasio Potasio 0.0050

2.5 Interpretación

Biomasa (gr). Media 1082.17 g y mediana 991.83 g. La media supera a la mediana y la asimetría es positiva (0.94), lo que indica una distribución sesgada a la derecha: la mayoría de las muestras producen biomasas moderadas y unas pocas alcanzan valores muy altos (máximo 2337.33 g frente a un mínimo de 369.82 g). Es la variable más dispersa en términos relativos (CV = 50.5 %), lo cual es esperable porque proviene de ambientes muy distintos. El boxplot no señala atípicos (el rango intercuartílico es amplio, 692.06 g), pero la prueba de Shapiro-Wilk rechaza la normalidad (p = 0.0008), congruente con el sesgo observado.

pH. Varía entre 3.20 y 7.45, con media 4.61 y mediana 4.45. Se trata de suelos predominantemente ácidos: el 75 % de las muestras tiene pH ≤ 5.35. La asimetría positiva (0.90) se explica por un pequeño grupo de muestras con pH neutro a ligeramente alcalino (≈ 7.1–7.45), que se separa del resto. El histograma sugiere una distribución bimodal, no una sola población homogénea.

Sal_suelo (salinidad del suelo). Es la variable más homogénea de todas (CV = 12.3 %, DE = 3.72). Media 30.27 y mediana 30.00 prácticamente coinciden y la asimetría es baja (0.32), de modo que su distribución es aproximadamente simétrica. El rango es estrecho (24 a 38) y Shapiro-Wilk no rechaza la normalidad al 5 % (p = 0.0507, en el límite). Esta poca variabilidad anticipa que difícilmente explicará las diferencias de biomasa.

Zinc. Media 17.83 y mediana 19.24, con asimetría negativa (−0.69). El boxplot marca cinco valores atípicos bajos (0.21, 0.26, 0.27, 0.32 y 0.37), muy alejados del grueso de los datos (Q1 = 13.99). No son errores de registro: corresponden justamente a las muestras de pH más alto, y es un hecho conocido en química de suelos que la disponibilidad del zinc cae drásticamente cuando el pH sube hacia la neutralidad. Por eso se conservan en el análisis, pero se reportan explícitamente.

Potasio. Media 797.38 y mediana 773.30, ligeramente asimétrica a la derecha (0.50), con amplia dispersión (350.73 a 1441.67; CV = 37.3 %). Sin valores atípicos según el criterio de 1.5·RIC.

Síntesis: ninguna variable, salvo posiblemente Sal_suelo, se comporta como normal; varias muestran sesgo y la estructura bimodal del pH sugiere que las muestras provienen de ambientes contrastantes. Esto justifica pasar al análisis bivariado.

3 Punto (b) — Análisis exploratorio bivariado

El objetivo es: cuantificar y visualizar la relación entre la biomasa y cada covariable, para identificar cuál se asocia más fuertemente con la respuesta.

Se usa el coeficiente de correlación de Pearson (r), que mide la fuerza y el sentido de la relación lineal: su signo indica la dirección y su valor absoluto, entre 0 y 1, la intensidad. Se acompaña con el coeficiente de Spearman, que detecta relaciones monótonas aunque no sean lineales, y con los diagramas de dispersión, que son los que realmente permiten juzgar la forma de la relación.

covariables <- c("pH", "Sal_suelo", "Zinc", "Potasio")

tabla_cor <- data.frame(
  Covariable = covariables,
  r_Pearson  = sapply(covariables, function(v) cor(Salinidad$Biomasa, Salinidad[[v]])),
  p_valor    = sapply(covariables, function(v) cor.test(Salinidad$Biomasa, Salinidad[[v]])$p.value),
  r2         = sapply(covariables, function(v) cor(Salinidad$Biomasa, Salinidad[[v]])^2),
  rho_Spearman = sapply(covariables, function(v)
                   cor(Salinidad$Biomasa, Salinidad[[v]], method = "spearman")),
  row.names = NULL
)
kable(tabla_cor, digits = 4,
      caption = "Correlación entre la biomasa y cada covariable del suelo")
Correlación entre la biomasa y cada covariable del suelo
Covariable r_Pearson p_valor r2 rho_Spearman
pH 0.9281 0.0000 0.8614 0.8777
Sal_suelo -0.0666 0.6639 0.0044 -0.1571
Zinc -0.7815 0.0000 0.6107 -0.6250
Potasio -0.0732 0.6328 0.0054 0.1016 
kable(round(cor(Salinidad), 3), caption = "Matriz de correlación entre todas las variables")
Matriz de correlación entre todas las variables
Biomasa pH Sal_suelo Zinc Potasio
Biomasa 1.000 0.928 -0.067 -0.781 -0.073
pH 0.928 1.000 -0.045 -0.720 0.032
Sal_suelo -0.067 -0.045 1.000 -0.427 -0.020
Zinc -0.781 -0.720 -0.427 1.000 0.079
Potasio -0.073 0.032 -0.020 0.079 1.000 
Salinidad %>%
  pivot_longer(all_of(covariables), names_to = "Covariable", values_to = "Valor") %>%
  ggplot(aes(x = Valor, y = Biomasa)) +
  geom_point(color = "steelblue", size = 2, alpha = 0.8) +
  geom_smooth(method = "lm", se = TRUE, color = "darkred", linewidth = 0.8) +
  facet_wrap(~ Covariable, scales = "free_x", ncol = 2) +
  labs(title = "Biomasa frente a cada covariable del suelo",
       x = NULL, y = "Biomasa (gr)") +
  theme_bw()

3.1 Interpretación

pH — r = 0.928 (p < 0.001). Correlación positiva y muy fuerte. El diagrama de dispersión muestra una nube claramente ascendente y bastante ajustada a la recta: a mayor pH del suelo, mayor producción de biomasa. Su r² = 0.861 indica que el pH explica por sí solo cerca del 86 % de la variabilidad de la biomasa. Es, con diferencia, la variable más asociada con la respuesta.

Zinc — r = −0.782 (p < 0.001). Correlación negativa y fuerte (r² = 0.611): las muestras con más zinc disponible tienden a producir menos biomasa. Conviene interpretarla con cautela, porque el zinc está a su vez fuertemente correlacionado con el pH (r = −0.72 entre ambas), de modo que parte de este efecto es un reflejo indirecto del pH y no necesariamente un efecto propio del zinc.

Sal_suelo — r = −0.067 (p = 0.664). Correlación prácticamente nula y no significativa. La nube de puntos no muestra ningún patrón. Coherente con lo visto en (a): su rango de variación es demasiado estrecho para generar diferencias.

Potasio — r = −0.073 (p = 0.633). También sin relación apreciable con la biomasa.

Conclusión: las variables con mayor relación con la biomasa son el pH (la más fuerte, positiva) y el Zinc (fuerte, negativa). La Sal_suelo y el Potasio no presentan asociación significativa. La variable seleccionada para el literal (c) es el pH, por tener la correlación de mayor magnitud.

4 Punto (c) — ANOVA de una vía sobre el pH categorizado

4.1 Categorización en tres niveles

El pH es una variable continua, mientras que un ANOVA de una vía requiere un factor categórico. Por eso se discretiza en tres niveles.

Criterio elegido: terciles. Se cortan los datos en los percentiles 33.3 y 66.7, lo que produce tres grupos de tamaño 15 cada uno. La justificación es:

  1. Es un criterio empírico y no arbitrario: nace de la distribución observada y no de puntos de corte impuestos desde fuera.
  2. Genera un diseño perfectamente balanceado (n = 15 por nivel), lo cual maximiza la potencia del ANOVA y lo vuelve más robusto frente a posibles desviaciones de los supuestos.
Cortes de terciles usados para categorizar el pH
0% 33.33333% 66.66667% 100%
3.2 3.883 4.9 7.45
Tamaño y rango de pH por nivel
pH_cat n pH_min pH_max
Bajo 15 3.20 3.75
Medio 15 3.95 4.85
Alto 15 5.00 7.45

Los niveles quedan definidos así: Bajo = pH 3.20–3.75, Medio = pH 3.95–4.85, Alto = pH 5.00–7.45.

4.2 Estadísticas descriptivas por nivel

Salinidad %>%
  group_by(pH_cat) %>%
  summarise(n = n(),
            Media = mean(Biomasa), DE = sd(Biomasa),
            CV_pct = 100 * sd(Biomasa) / mean(Biomasa),
            Minimo = min(Biomasa), Maximo = max(Biomasa)) %>%
  kable(digits = 2, caption = "Biomasa (gr) según el nivel de pH")
Biomasa (gr) según el nivel de pH
pH_cat n Media DE CV_pct Minimo Maximo
Bajo 15 593.02 164.83 27.79 369.82 977.52
Medio 15 1048.11 244.91 23.37 568.46 1491.28
Alto 15 1605.39 547.60 34.11 765.28 2337.33 
ggplot(Salinidad, aes(x = pH_cat, y = Biomasa, fill = pH_cat)) +
  geom_boxplot(alpha = 0.7, outlier.color = "red") +
  geom_jitter(width = 0.12, alpha = 0.6, size = 1.8) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "black") +
  labs(title = "Biomasa según el nivel de pH del suelo",
       subtitle = "El rombo negro indica la media de cada grupo",
       x = "Nivel de pH", y = "Biomasa (gr)") +
  scale_fill_brewer(palette = "Blues") +
  theme_bw() + theme(legend.position = "none")

Las medias ascienden de forma marcada: 593.02 g (Bajo), 1048.11 g (Medio) y 1605.39 g (Alto). Nótese, sin embargo, que la desviación estándar también crece con la media (164.83 → 244.91 → 547.60). Esto es una señal temprana de posible heterogeneidad de varianzas, que se verificará formalmente.

4.3 Modelo ANOVA sobre la variable original

El modelo es:

\[ Y_{ij} = \mu + \tau_i + \varepsilon_{ij}, \qquad i = 1,2,3 \ (\text{niveles}); \quad j = 1,\dots,15 \]

donde \(Y_{ij}\) es la biomasa de la muestra \(j\) del nivel \(i\), \(\mu\) la media general, \(\tau_i\) el efecto del nivel \(i\) de pH y \(\varepsilon_{ij} \sim N(0,\sigma^2)\) independientes.

Las hipótesis son:

  • \(H_0: \tau_1 = \tau_2 = \tau_3 = 0\) (el nivel de pH no afecta la biomasa media)
  • \(H_1:\) al menos un \(\tau_i \neq 0\)
modelo0 <- aov(Biomasa ~ pH_cat, data = Salinidad)

tabla_anova0 <- as.data.frame(summary(modelo0)[[1]])
tabla_anova0 <- cbind(Fuente = trimws(rownames(tabla_anova0)), tabla_anova0)
rownames(tabla_anova0) <- NULL
kable(tabla_anova0, digits = c(0, 0, 2, 2, 2, 10),
      col.names = c("Fuente", "gl", "SC", "CM", "F", "p-valor"),
      caption = "ANOVA de una vía — Biomasa ~ nivel de pH (escala original)")
ANOVA de una vía — Biomasa ~ nivel de pH (escala original)
Fuente gl SC CM F p-valor
pH_cat 2 7712683 3856341.6 29.89 8.4e-09
Residuals 42 5418235 129005.6 NA NA

El estadístico es F = 29.89 con 2 y 42 grados de libertad, p ≈ 8.45 × 10⁻⁹, así que el efecto es altamente significativo. Pero antes de darlo por válido hay que verificar los supuestos.

4.4 Verificación de supuestos (modelo original)

Los supuestos se evalúan sobre los residuales, no sobre los datos crudos.

par(mfrow = c(1, 2))
plot(modelo0, which = 2, main = "QQ-plot de residuales")
plot(modelo0, which = 1, main = "Residuales vs ajustados")

par(mfrow = c(1, 1))
sh0 <- shapiro.test(residuals(modelo0))
lv0 <- leveneTest(Biomasa ~ pH_cat, data = Salinidad)
bt0 <- bartlett.test(Biomasa ~ pH_cat, data = Salinidad)

tabla_supuestos0 <- data.frame(
  Prueba = c("Shapiro-Wilk (normalidad)", "Levene (varianzas)", "Bartlett (varianzas)"),
  Estadistico = c(sh0$statistic, lv0$`F value`[1], bt0$statistic),
  p_valor = c(sh0$p.value, lv0$`Pr(>F)`[1], bt0$p.value)
)
kable(tabla_supuestos0, digits = 4,
      caption = "Verificación de supuestos — modelo sobre Biomasa original")
Verificación de supuestos — modelo sobre Biomasa original
Prueba Estadistico p_valor
W Shapiro-Wilk (normalidad) 0.9732 0.3749
Levene (varianzas) 8.6753 0.0007
Bartlett’s K-squared Bartlett (varianzas) 20.0843 0.0000
  • Normalidad: se cumple. Shapiro-Wilk sobre los residuales da p = 0.3749 > 0.05, por lo que no se rechaza \(H_0\); el QQ-plot muestra los puntos alineados sobre la diagonal.
  • Homogeneidad de varianzas: NO se cumple. La prueba de Levene da p = 0.0007 < 0.05 (Bartlett, aún más drástica, p < 0.001), de modo que se rechaza la igualdad de varianzas. El gráfico de residuales contra valores ajustados lo confirma: la dispersión se abre en forma de embudo a medida que crece la media.

Decisión: el incumplimiento de homocedasticidad invalida el uso directo de la prueba F y, sobre todo, de la LSD (que usa un único cuadrado medio del error común a todos los grupos). Como la varianza crece proporcionalmente con la media —patrón típico de datos de biomasa—, se aplica la transformación logarítmica, que es el remedio estándar en esta situación.

4.5 Modelo corregido: ANOVA sobre log(Biomasa)

Salinidad$logBiomasa <- log(Salinidad$Biomasa)

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

tabla_anova <- as.data.frame(summary(modelo)[[1]])
tabla_anova <- cbind(Fuente = trimws(rownames(tabla_anova)), tabla_anova)
rownames(tabla_anova) <- NULL
kable(tabla_anova, digits = c(0, 0, 4, 4, 2, 10),
      col.names = c("Fuente", "gl", "SC", "CM", "F", "p-valor"),
      caption = "ANOVA de una vía — log(Biomasa) ~ nivel de pH")
ANOVA de una vía — log(Biomasa) ~ nivel de pH
Fuente gl SC CM F p-valor
pH_cat 2 7.1536 3.5768 40.81 1e-10
Residuals 42 3.6807 0.0876 NA NA 
par(mfrow = c(1, 2))
plot(modelo, which = 2, main = "QQ-plot (modelo log)")
plot(modelo, which = 1, main = "Residuales vs ajustados (modelo log)")

par(mfrow = c(1, 1))
sh1 <- shapiro.test(residuals(modelo))
lv1 <- leveneTest(logBiomasa ~ pH_cat, data = Salinidad)

tabla_supuestos1 <- data.frame(
  Prueba = c("Shapiro-Wilk (normalidad)", "Levene (varianzas)"),
  Estadistico = c(sh1$statistic, lv1$`F value`[1]),
  p_valor = c(sh1$p.value, lv1$`Pr(>F)`[1])
)
kable(tabla_supuestos1, digits = 4,
      caption = "Verificación de supuestos — modelo sobre log(Biomasa)")
Verificación de supuestos — modelo sobre log(Biomasa)
Prueba Estadistico p_valor
W Shapiro-Wilk (normalidad) 0.9790 0.5780
Levene (varianzas) 1.4135 0.2546

Con la transformación se cumplen ambos supuestos:

  • Normalidad de residuales: Shapiro-Wilk p = 0.578 > 0.05. ✔
  • Homogeneidad de varianzas: Levene p = 0.2546 > 0.05. ✔

Y el resultado del ANOVA es F = 40.81 (gl = 2, 42), p ≈ 1.42 × 10⁻¹⁰. Se rechaza \(H_0\): el nivel de pH del suelo genera diferencias altamente significativas en la biomasa producida. Nótese que la conclusión coincide con la del modelo sin transformar, lo que refuerza su solidez; la diferencia es que ahora la inferencia es estadísticamente válida.

4.6 Comparaciones post-hoc: prueba de LSD

El ANOVA solo indica que al menos un nivel difiere; no dice cuáles. La prueba LSD (Diferencia Mínima Significativa) de Fisher compara todos los pares calculando la diferencia mínima que debe existir entre dos medias para considerarlas distintas:

\[ \text{LSD} = t_{\alpha/2,\,gl_{error}} \sqrt{CME\left(\frac{1}{n_i}+\frac{1}{n_j}\right)} \]

Con \(CME = 0.0876\), \(gl_{error} = 42\), \(t_{0.025,42} = 2.018\) y \(n_i = n_j = 15\):

\[ \text{LSD} = 2.018 \sqrt{0.0876\left(\tfrac{1}{15}+\tfrac{1}{15}\right)} = 0.2181 \]

Es decir, dos niveles se declaran diferentes si sus medias en escala logarítmica difieren en más de 0.2181.

lsd <- LSD.test(modelo, "pH_cat", alpha = 0.05, console = FALSE)

plot(lsd, main = "Prueba LSD: agrupación de los niveles de pH",
     xlab = "Nivel de pH", ylab = "log(Biomasa)")

# Medias en escala log y su retransformación a gramos (medias geométricas)
medias_log <- tapply(Salinidad$logBiomasa, Salinidad$pH_cat, mean)
kable(data.frame(Nivel = names(medias_log),
           Media_log = round(as.numeric(medias_log), 3),
           Media_geometrica_gr = round(exp(as.numeric(medias_log)), 1)),
      caption = "Medias por nivel (escala log y retransformadas a gramos)")
Medias por nivel (escala log y retransformadas a gramos)
Nivel Media_log Media_geometrica_gr
Bajo 6.352 573.5
Medio 6.927 1019.4
Alto 7.323 1514.5

4.6.1 Interpretación de la prueba LSD

Comparación Diferencia (log) LSD p-valor ¿Difieren?
Alto vs Bajo 0.9711 0.2181 2.49×10⁻¹¹ Sí
Alto vs Medio 0.3959 0.2181 0.00069 Sí
Medio vs Bajo 0.5752 0.2181 3.73×10⁻⁶ Sí

Los tres niveles difieren significativamente entre sí, ya que las tres diferencias superan el valor LSD de 0.2181 y los tres p-valores son menores a 0.05. En consecuencia, la prueba asigna una letra distinta a cada nivel (a = Alto, b = Medio, c = Bajo): no hay ningún par de grupos que pueda considerarse equivalente.

Retransformando a la escala original, las medias geométricas de biomasa son aproximadamente 573.6 g (pH Bajo), 1019.4 g (pH Medio) y 1514.7 g (pH Alto). El paso de pH bajo a medio casi duplica la producción, y el paso a pH alto la multiplica por cerca de 2.6 respecto al nivel más ácido.

5 Conclusiones del Punto 1

  1. El análisis univariado mostró variables con comportamientos muy distintos: la Biomasa y el pH son asimétricas a la derecha, el Zinc presenta cinco valores atípicos bajos asociados a los suelos de pH más alto, y la Sal_suelo es la variable más homogénea de todas (CV = 12.3 %).

  2. El análisis bivariado identificó al pH como la covariable más asociada con la biomasa (r = 0.928, r² = 0.861), seguida del Zinc en sentido inverso (r = −0.782). La Sal_suelo (r = −0.067) y el Potasio (r = −0.073) no mostraron relación significativa.

  3. Categorizando el pH en terciles, el ANOVA de una vía detectó diferencias altamente significativas en la biomasa (F = 40.81; p ≈ 1.42 × 10⁻¹⁰). Fue necesario transformar la respuesta con logaritmo, porque el modelo sobre la escala original incumplía la homogeneidad de varianzas (Levene p = 0.0007); tras la transformación ambos supuestos se satisfacen.

  4. La prueba LSD indicó que los tres niveles de pH difieren entre sí, con un incremento monótono de la biomasa a medida que el suelo se acerca a la neutralidad. Desde el punto de vista agronómico, esto sugiere que la acidez del suelo es el principal factor limitante de la producción de esta forrajera, y que prácticas de encalado orientadas a elevar el pH tendrían un impacto directo sobre el rendimiento.

  5. Advertencia sobre la causalidad: se trata de un estudio observacional, no de un experimento con asignación aleatoria de tratamientos. El pH está correlacionado con otras propiedades del suelo (notablemente el Zinc, r = −0.72), de modo que el efecto estimado puede estar parcialmente confundido con ellas. Las diferencias detectadas son reales, pero atribuirlas exclusivamente al pH requeriría un ensayo controlado.