Para estudiar la relación entre características del suelo y la producción de biomasa de una planta forrajera natural se obtuvieron 45 muestras en diferentes ambientes. En cada una se midió la biomasa (variable respuesta, en gramos) y cuatro covariables del suelo: pH, Salinidad, Zinc y Potasio.
# Ajusta la ruta según donde tengas el archivo
load("C:/Users/isapa/Downloads/Salinidad.RData")
df <- Salinidad
str(df)
## 'data.frame': 45 obs. of 5 variables:
## $ Biomasa : num 765 954 828 755 896 ...
## $ pH : num 5 4.7 4.2 4.4 5.55 5.5 4.25 4.45 4.75 4.6 ...
## $ Salinidad: int 33 35 32 30 33 33 36 30 38 30 ...
## $ Zinc : num 16.5 14 15.3 17.3 22.3 ...
## $ Potasio : num 1442 1299 1154 1045 522 ...
summary(df)
## Biomasa pH Salinidad Zinc
## Min. : 369.8 Min. :3.200 Min. :24.00 Min. : 0.2105
## 1st Qu.: 654.8 1st Qu.:3.450 1st Qu.:27.00 1st Qu.:13.9852
## Median : 991.8 Median :4.450 Median :30.00 Median :19.2420
## Mean :1082.2 Mean :4.609 Mean :30.27 Mean :17.8308
## 3rd Qu.:1346.9 3rd Qu.:5.350 3rd Qu.:33.00 3rd Qu.:22.6758
## Max. :2337.3 Max. :7.450 Max. :38.00 Max. :31.2865
## Potasio
## Min. : 350.7
## 1st Qu.: 527.0
## Median : 773.3
## Mean : 797.4
## 3rd Qu.: 954.1
## Max. :1441.7
descriptivos <- sapply(df, function(x) c(
media = mean(x), de = sd(x), cv = sd(x) / mean(x) * 100,
min = min(x), max = max(x)
))
round(descriptivos, 2)
## Biomasa pH Salinidad Zinc Potasio
## media 1082.17 4.61 30.27 17.83 797.38
## de 546.29 1.25 3.72 8.27 297.58
## cv 50.48 27.22 12.29 46.40 37.32
## min 369.82 3.20 24.00 0.21 350.73
## max 2337.33 7.45 38.00 31.29 1441.67
skew <- function(x) { n <- length(x); m <- mean(x); s <- sd(x); (sum((x - m)^3) / n) / s^3 }
kurt <- function(x) { n <- length(x); m <- mean(x); s <- sd(x); (sum((x - m)^4) / n) / s^4 - 3 }
round(sapply(df, function(x) c(asimetria = skew(x), curtosis = kurt(x))), 2)
## Biomasa pH Salinidad Zinc Potasio
## asimetria 0.91 0.87 0.31 -0.66 0.48
## curtosis -0.02 -0.03 -1.04 -0.06 -1.00
for (v in names(df)) {
b <- boxplot.stats(df[[v]])
cat(v, "-> outliers:",
if (length(b$out) > 0) paste(round(b$out, 2), collapse = ", ") else "ninguno", "\n")
}
## Biomasa -> outliers: ninguno
## pH -> outliers: ninguno
## Salinidad -> outliers: ninguno
## Zinc -> outliers: 0.37, 0.27, 0.32, 0.26, 0.21
## Potasio -> outliers: ninguno
for (v in names(df)) {
p <- ggplot(df, aes(x = .data[[v]])) +
geom_histogram(bins = 10, fill = "#4C72B0", color = "white") +
labs(title = paste("Histograma -", v), x = v, y = "Frecuencia") +
theme_minimal(base_size = 13)
print(p)
}
for (v in names(df)) {
p <- ggplot(df, aes(y = .data[[v]])) +
geom_boxplot(fill = "#DD8452") +
labs(title = paste("Boxplot -", v), y = v) +
theme_minimal(base_size = 13) +
theme(axis.text.x = element_blank())
print(p)
}
Las cinco variables analizadas muestran comportamientos muy distintos entre sí, tanto en su forma de distribución como en su nivel de dispersión relativa. Biomasa, pH y Potasio presentan sesgo positivo (asimetrías de 0.91, 0.87 y 0.48 respectivamente), lo que indica que en las tres variables la mayoría de las observaciones se concentra en valores bajos o moderados, mientras que un subgrupo más pequeño de muestras presenta valores notablemente más altos que estiran la cola de la distribución hacia la derecha. En el caso de Biomasa, esto se traduce en que la mayoría de los 45 sitios produjeron entre 400 y 1500 g, pero un grupo reducido de sitios alcanzó valores entre 1900 y 2337 g — un patrón que, como se confirma en el literal (b), está directamente asociado a que esos mismos sitios presentan también los valores más altos de pH.
Salinidad se comporta de forma distinta: es la variable más simétrica del conjunto (asimetría de 0.31, la más cercana a cero) y además la de menor dispersión relativa (CV = 12.3%, frente a valores entre 27% y 50% en las demás variables). Esto sugiere que, a diferencia de las otras covariables, la salinidad del suelo varió relativamente poco entre los 45 ambientes muestreados — un indicio temprano de que difícilmente iba a mostrar una relación fuerte con la biomasa, algo que efectivamente se confirmó en el análisis bivariado.
Zinc es la única variable con sesgo negativo (asimetría de -0.66): aquí la mayoría de las muestras presenta niveles moderados a altos (10-31), pero cinco observaciones muestran valores casi nulos (0.21-0.37), identificadas como atípicas mediante el criterio de rango intercuartílico. Este subgrupo diferenciado podría reflejar un tipo particular de suelo con deficiencia real de zinc, o bien una limitación del método de medición (valores por debajo del límite de detección del instrumento); en cualquier caso, no se eliminaron del análisis porque no hay evidencia de que sean errores de captura, y descartarlos sin justificación reduciría artificialmente la variabilidad real del suelo que el estudio busca capturar.
En conjunto, este análisis univariado no es solo un paso descriptivo previo: la asimetría compartida entre Biomasa y pH, y el subgrupo diferenciado en Zinc, ya anticipan patrones que se confirman estadísticamente en los literales siguientes.
covariables <- c("pH", "Salinidad", "Zinc")
cat("--- Pearson ---\n")
## --- Pearson ---
for (v in covariables) {
ct <- cor.test(df$Biomasa, df[[v]], method = "pearson")
cat(sprintf("%-10s r = %.3f p = %.4f\n", v, ct$estimate, ct$p.value))
}
## pH r = 0.928 p = 0.0000
## Salinidad r = -0.067 p = 0.6639
## Zinc r = -0.781 p = 0.0000
cat("\n--- Spearman ---\n")
##
## --- Spearman ---
for (v in covariables) {
ct <- suppressWarnings(cor.test(df$Biomasa, df[[v]], method = "spearman"))
cat(sprintf("%-10s rho = %.3f p = %.4f\n", v, ct$estimate, ct$p.value))
}
## pH rho = 0.878 p = 0.0000
## Salinidad rho = -0.157 p = 0.3027
## Zinc rho = -0.625 p = 0.0000
for (v in covariables) {
p <- ggplot(df, aes(x = .data[[v]], y = Biomasa)) +
geom_point(color = "#4C72B0", size = 2.2, alpha = 0.8) +
geom_smooth(method = "lm", se = TRUE, color = "#C44E52") +
labs(title = paste("Biomasa vs", v)) +
theme_minimal(base_size = 13)
print(p)
}
Para identificar qué covariable se relaciona más con la biomasa, se calcularon coeficientes de correlación de Pearson y de Spearman entre Biomasa y cada una de las tres covariables (pH, Salinidad y Zinc). Se utilizaron ambos coeficientes de manera deliberada y no solo Pearson, precisamente porque el análisis univariado había evidenciado sesgo en Biomasa y en Zinc: Pearson, al calcularse sobre los valores originales, es sensible a esa asimetría y a valores atípicos, mientras que Spearman, al trabajar sobre los rangos de los datos en lugar de sus valores absolutos, resulta más robusto ante ese tipo de distorsión. Que ambos coeficientes hayan arrojado conclusiones consistentes —mismo orden de fuerza y mismo signo en las tres covariables— es evidencia de que las relaciones encontradas son reales y no producto de un puñado de valores extremos.
Los resultados muestran a pH como la covariable con mayor relación con la biomasa (r = 0.928, ρ = 0.878), una correlación positiva y notablemente fuerte: a mayor pH del suelo, mayor producción de biomasa, y el diagrama de dispersión confirma que este patrón sigue de forma consistente una tendencia lineal a lo largo de todo el rango de pH observado, sin desviaciones marcadas de ese comportamiento. Zinc presenta la segunda relación más fuerte, pero de signo negativo (r = -0.781, ρ = -0.625): a mayor concentración de zinc en el suelo, menor biomasa producida. Salinidad, en cambio, no muestra relación relevante con la biomasa (r = -0.067, ρ = -0.157, ambos no significativos), lo cual es coherente con lo observado en el análisis univariado, donde esta variable presentó la menor variabilidad relativa del conjunto: al variar poco entre los distintos ambientes muestreados, es esperable que su capacidad de explicar diferencias en biomasa sea limitada.
Con base en estos resultados, se seleccionó pH como la covariable con mayor relación con la biomasa para desarrollar el análisis de varianza del literal (c), al ser la única de las tres cuya correlación es simultáneamente la más fuerte en magnitud y la más consistente entre ambos coeficientes.
terciles <- quantile(df$pH, probs = c(1/3, 2/3))
terciles
## 33.33333% 66.66667%
## 3.883333 4.900000
df$pH_nivel <- cut(df$pH,
breaks = c(-Inf, terciles[1], terciles[2], Inf),
labels = c("Bajo", "Medio", "Alto"))
df$pH_nivel <- factor(df$pH_nivel, levels = c("Bajo", "Medio", "Alto"))
table(df$pH_nivel)
##
## Bajo Medio Alto
## 15 15 15
aggregate(pH ~ pH_nivel, df, function(x) c(min = min(x), max = max(x)))
## pH_nivel pH.min pH.max
## 1 Bajo 3.20 3.75
## 2 Medio 3.95 4.85
## 3 Alto 5.00 7.45
aggregate(Biomasa ~ pH_nivel, df, function(x) c(media = mean(x), de = sd(x)))
## pH_nivel Biomasa.media Biomasa.de
## 1 Bajo 593.0233 164.8299
## 2 Medio 1048.1093 244.9093
## 3 Alto 1605.3853 547.6014
ggplot(df, aes(x = pH_nivel, y = Biomasa, fill = pH_nivel)) +
geom_boxplot(alpha = 0.85, show.legend = FALSE) +
geom_jitter(width = 0.08, alpha = 0.5, size = 1.8) +
scale_fill_manual(values = c("#8FBCDB", "#4C72B0", "#1B3B6F")) +
labs(title = "Biomasa según nivel de pH (terciles)",
x = "Nivel de pH", y = "Biomasa (g)") +
theme_minimal(base_size = 13)
mod <- aov(Biomasa ~ pH_nivel, data = df)
summary(mod)
## Df Sum Sq Mean Sq F value Pr(>F)
## pH_nivel 2 7712683 3856342 29.89 8.45e-09 ***
## Residuals 42 5418235 129006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res <- residuals(mod)
shapiro.test(res)
##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.97316, p-value = 0.3749
qqnorm(res, main = "QQ-plot de residuales")
qqline(res, col = "red")
leveneTest(Biomasa ~ pH_nivel, data = df)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 8.6753 0.000702 ***
## 42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bartlett.test(Biomasa ~ pH_nivel, data = df)
##
## Bartlett test of homogeneity of variances
##
## data: Biomasa by pH_nivel
## Bartlett's K-squared = 20.084, df = 2, p-value = 4.353e-05
La covariable pH se categorizó en tres niveles (Bajo, Medio, Alto) usando terciles, de modo que cada grupo quedara compuesto por 15 observaciones, garantizando grupos balanceados para el análisis de varianza. Las medias de biomasa resultantes muestran un incremento claro y consistente entre niveles: 593 g en el nivel Bajo, 1048 g en el nivel Medio y 1605 g en el nivel Alto — un patrón que replica, ahora de forma categórica, la relación lineal positiva ya identificada en el literal (b).
El ANOVA de una vía confirmó que estas diferencias son estadísticamente significativas (F = 29.89, p < 0.0001): la variación de biomasa explicada por el nivel de pH es considerablemente mayor que la variación residual dentro de cada grupo, por lo que se rechaza la hipótesis de que las tres medias son iguales. Sin embargo, la validez de esta prueba depende de que se cumplan sus supuestos, por lo que estos se verificaron explícitamente antes de interpretar el resultado como concluyente. La normalidad de los residuales del modelo (no de la variable Biomasa cruda, sino de los residuales, que son los que efectivamente exige el supuesto) se evaluó con la prueba de Shapiro-Wilk, sin evidencia para rechazarla (p = 0.375), lo cual fue corroborado visualmente mediante el gráfico QQ, donde los puntos siguen de forma razonable la línea de normalidad teórica. En contraste, el supuesto de homogeneidad de varianzas no se cumplió: tanto la prueba de Levene (p = 0.0007) como la de Bartlett (p < 0.0001) rechazaron la homocedasticidad, resultado consistente con las desviaciones estándar observadas por grupo (165 en el nivel Bajo, 245 en el Medio y 548 en el Alto), donde la dispersión de la biomasa crece claramente a medida que aumenta el nivel de pH.
Esta violación del supuesto es relevante porque tanto el ANOVA clásico como la prueba post-hoc LSD solicitada en el enunciado calculan un error estándar común a partir de una varianza combinada de los tres grupos (el cuadrado medio del error), lo cual deja de ser representativo cuando las varianzas de los grupos son tan distintas entre sí. Por esta razón, además de aplicar la prueba LSD como lo pedía el ejercicio, se contrastaron sus resultados con una alternativa robusta a la heterocedasticidad (ANOVA de Welch y prueba post-hoc de Games-Howell), que no asume varianzas iguales entre grupos.
s <- summary(mod)[[1]]
MSE <- s["Residuals", "Mean Sq"]
dfres <- s["Residuals", "Df"]
alpha <- 0.05
tcrit <- qt(1 - alpha / 2, dfres)
medias <- aggregate(Biomasa ~ pH_nivel, df, mean)
ns <- table(df$pH_nivel)
niveles <- as.character(medias$pH_nivel)
pares <- combn(niveles, 2)
for (i in 1:ncol(pares)) {
g1 <- pares[1, i]; g2 <- pares[2, i]
m1 <- medias$Biomasa[medias$pH_nivel == g1]
m2 <- medias$Biomasa[medias$pH_nivel == g2]
n1 <- ns[g1]; n2 <- ns[g2]
diff <- m1 - m2
se <- sqrt(MSE * (1/n1 + 1/n2))
lsd <- tcrit * se
tstat <- diff / se
pval <- 2 * (1 - pt(abs(tstat), dfres))
cat(sprintf("%s - %s: diff = %.1f | LSD = %.1f | t = %.2f | p = %.5f | %s\n",
g1, g2, diff, lsd, tstat, pval,
ifelse(abs(diff) > lsd, "DIFIEREN", "no difieren")))
}
## Bajo - Medio: diff = -455.1 | LSD = 264.7 | t = -3.47 | p = 0.00122 | DIFIEREN
## Bajo - Alto: diff = -1012.4 | LSD = 264.7 | t = -7.72 | p = 0.00000 | DIFIEREN
## Medio - Alto: diff = -557.3 | LSD = 264.7 | t = -4.25 | p = 0.00012 | DIFIEREN
oneway.test(Biomasa ~ pH_nivel, data = df, var.equal = FALSE)
##
## One-way analysis of means (not assuming equal variances)
##
## data: Biomasa and pH_nivel
## F = 34.727, num df = 2.000, denom df = 24.736, p-value = 6.58e-08
games_howell <- function(x1, x2, n1, n2) {
m1 <- mean(x1); m2 <- mean(x2); v1 <- var(x1); v2 <- var(x2)
se <- sqrt(v1 / n1 + v2 / n2)
tstat <- (m1 - m2) / se
df_gh <- (v1/n1 + v2/n2)^2 / ((v1/n1)^2/(n1-1) + (v2/n2)^2/(n2-1))
p <- ptukey(sqrt(2) * abs(tstat), 3, df_gh, lower.tail = FALSE)
c(diff = m1 - m2, t = tstat, df = df_gh, p = p)
}
for (i in 1:2) for (j in (i + 1):3) {
g1 <- niveles[i]; g2 <- niveles[j]
x1 <- df$Biomasa[df$pH_nivel == g1]
x2 <- df$Biomasa[df$pH_nivel == g2]
r <- games_howell(x1, x2, length(x1), length(x2))
cat(sprintf("%s - %s: diff = %.1f t = %.2f df = %.1f p = %.5f\n",
g1, g2, r["diff"], r["t"], r["df"], r["p"]))
}
## Bajo - Medio: diff = -455.1 t = -5.97 df = 24.5 p = 0.00001
## Bajo - Alto: diff = -1012.4 t = -6.86 df = 16.5 p = 0.00001
## Medio - Alto: diff = -557.3 t = -3.60 df = 19.4 p = 0.00507
La prueba LSD indicó diferencias significativas entre los tres pares de niveles (Bajo–Medio: diferencia de 455 g, p = 0.0012; Bajo–Alto: 1012 g, p < 0.0001; Medio–Alto: 557 g, p = 0.0001), y estas mismas conclusiones se mantuvieron con Games-Howell (p = 0.00001, 0.00001 y 0.00507, respectivamente). Que ambos métodos —uno que asume varianzas iguales y otro que no— coincidan en que los tres niveles difieren significativamente entre sí, con el mismo orden de magnitud (Bajo < Medio < Alto), da respaldo a que el hallazgo es robusto y no depende de la violación del supuesto de homocedasticidad detectada.
En conclusión, existe evidencia estadística sólida de que el nivel de pH del suelo genera diferencias significativas en la producción de biomasa de la planta forrajera, con una relación de tipo creciente: a mayor nivel de pH, mayor biomasa, y los tres niveles categorizados son distinguibles entre sí. Cabe señalar, no obstante, una limitación metodológica: dado que pH mostró una relación lineal fuerte con la biomasa (r = 0.928), categorizarla en tres niveles discretos implica una simplificación que sacrifica parte de la información y de la potencia estadística disponible; un modelo de regresión lineal simple habría aprovechado de forma más completa la relación continua entre ambas variables. Se optó por el enfoque categórico porque era el solicitado en el ejercicio, pero es un punto que vale la pena reconocer como limitación del análisis y no como su forma óptima.