Efecto de la aplicación foliar de diferentes dosis de caolín (Al₂Si₂O₅(OH)₄) sobre el rendimiento en peso de los granos por planta de arveja (Pisum sativum)
Diseño en Bloques Completos al Azar (DBCA)
Razón de bloqueo: Humedad del suelo
July Acosta, Santiago Diaz, José Rincón, Paula Valbuena
1 Introducción
El presente documento muestra el análisis estadístico de un experimento establecido bajo un Diseño en Bloques Completos al Azar (DBCA), cuyo objetivo fue evaluar el efecto de diferentes tratamientos con caolín sobre el rendimiento de grano por planta de arveja (Pisum sativum).
La variable de respuesta evaluada fue el peso de grano por planta (g), considerando como factor de bloqueo la humedad del suelo con el fin de controlar parte de la variabilidad experimental.
1.1 Carga de librerías
#---------------------------------------------------------
# Instalación automática de paquetes
#---------------------------------------------------------
paquetes <- c(
"readxl",
"dplyr",
"ggplot2",
"agricolae",
"car",
"knitr"
)
instalar <- paquetes[!(paquetes %in% installed.packages()[,"Package"])]
if(length(instalar)>0){
install.packages(instalar)
}
lapply(paquetes, library, character.only = TRUE)## [[1]]
## [1] "readxl" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[2]]
## [1] "dplyr" "readxl" "stats" "graphics" "grDevices" "utils"
## [7] "datasets" "methods" "base"
##
## [[3]]
## [1] "ggplot2" "dplyr" "readxl" "stats" "graphics" "grDevices"
## [7] "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "agricolae" "ggplot2" "dplyr" "readxl" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[5]]
## [1] "car" "carData" "agricolae" "ggplot2" "dplyr" "readxl"
## [7] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [13] "base"
##
## [[6]]
## [1] "knitr" "car" "carData" "agricolae" "ggplot2" "dplyr"
## [7] "readxl" "stats" "graphics" "grDevices" "utils" "datasets"
## [13] "methods" "base"1.2 Bloqueos
Se observó que en el lote se presentó un gradiente de humedad bidireccional ocasionado por las zanjas que bordean el invernadero en los costados superior y derecho. Debido a esta variación ambiental, se establecieron bloques de acuerdo con el nivel de humedad del suelo, con el fin de reducir la variabilidad experimental y aumentar la precisión en la comparación de los tratamientos.
1.2.1 Gradiente de humedad
Figura 1. Gradiente bidireccional de humedad del área experimental.
1.2.2 Orden de tratamientos por bloque
Se organizaron los tratamientos de manera equitativa dentro de cada bloque, procurando que no coincidieran en el mismo surco ni en parcelas contiguas, con el fin de minimizar posibles efectos de vecindad y evitar sesgos en la evaluación de los tratamientos.
Figura 2. Organización final de los tratamientos por bloque
1.3 Importación de datos
## # A tibble: 6 × 4
## Humedad Planta Tratamiento Peso
## <chr> <dbl> <chr> <dbl>
## 1 Media alta 4 Control 31.3
## 2 Media alta 5 Control 31.9
## 3 Media alta 6 Control 6.34
## 4 Media alta 10 2,5% 28.6
## 5 Media alta 11 2,5% 54.6
## 6 Media alta 12 2,5% 34.2## tibble [36 × 4] (S3: tbl_df/tbl/data.frame)
## $ Humedad : chr [1:36] "Media alta" "Media alta" "Media alta" "Media alta" ...
## $ Planta : num [1:36] 4 5 6 10 11 12 16 17 18 7 ...
## $ Tratamiento: chr [1:36] "Control" "Control" "Control" "2,5%" ...
## $ Peso : num [1:36] 31.33 31.93 6.34 28.63 54.57 ...1.4 Preparación de los datos
1.5 Vista general
## Humedad Planta Tratamiento Peso
## Alta :9 Min. : 1.00 2,5% :12 Min. : 3.08
## Baja :9 1st Qu.: 9.75 5% :12 1st Qu.:13.03
## Media alta:9 Median :18.50 Control:12 Median :26.84
## Media baja:9 Mean :18.50 Mean :26.53
## 3rd Qu.:27.25 3rd Qu.:33.52
## Max. :36.00 Max. :77.071.6 Objetivo del análisis
Determinar si existen diferencias significativas en el rendimiento de grano por planta entre los tratamientos evaluados mediante un Diseño en Bloques Completos al Azar (DBCA), considerando la humedad del suelo como criterio de bloqueo.
Asimismo, verificar el cumplimiento de los supuestos del modelo y establecer cuáles tratamientos difieren estadísticamente mediante la prueba de comparaciones múltiples de Tukey.
2 Análisis descriptivo
El análisis descriptivo permite obtener una visión general del comportamiento de la variable respuesta (Peso) antes de realizar las pruebas de hipótesis.
2.1 Estadísticos descriptivos por tratamiento
#---------------------------------------------------------
# Estadísticos descriptivos
#---------------------------------------------------------
descriptivos <- datos %>%
group_by(Tratamiento) %>%
summarise(
n = n(),
Media = mean(Peso),
Mediana = median(Peso),
Desviacion = sd(Peso),
Minimo = min(Peso),
Maximo = max(Peso),
CV = (sd(Peso)/mean(Peso))*100
)
knitr::kable(
descriptivos,
digits = 2,
caption = "Estadísticos descriptivos por tratamiento"
)| Tratamiento | n | Media | Mediana | Desviacion | Minimo | Maximo | CV |
|---|---|---|---|---|---|---|---|
| 2,5% | 12 | 28.28 | 28.43 | 16.46 | 3.49 | 54.57 | 58.20 |
| 5% | 12 | 27.01 | 26.99 | 19.79 | 3.08 | 72.41 | 73.26 |
| Control | 12 | 24.30 | 19.65 | 18.68 | 6.34 | 77.07 | 76.86 |
2.2 Interpretación
cat("El tratamiento con mayor rendimiento promedio fue:",
descriptivos$Tratamiento[which.max(descriptivos$Media)],
"\n")## El tratamiento con mayor rendimiento promedio fue: 1cat("El tratamiento con menor rendimiento promedio fue:",
descriptivos$Tratamiento[which.min(descriptivos$Media)],
"\n")## El tratamiento con menor rendimiento promedio fue: 32.3 Distribución del rendimiento por tratamiento
ggplot(datos,
aes(x = Tratamiento,
y = Peso,
fill = Tratamiento))+
geom_boxplot(alpha=0.8)+
geom_jitter(width=0.15,
size=2,
colour="black")+
labs(
title="Distribución del rendimiento por tratamiento",
x="Tratamientos",
y="Peso de granos por planta (g)"
)+
theme_bw()+
theme(
legend.position="none",
plot.title=element_text(face="bold",hjust=0.5),
axis.title=element_text(face="bold")
)
2.4 Distribución del rendimiento por bloque
ggplot(datos,
aes(x = Humedad,
y = Peso,
fill = Humedad))+
geom_boxplot(alpha=0.8)+
geom_jitter(width=0.15,
size=2,
colour="black")+
labs(
title="Distribución del rendimiento por tratamiento",
x="Humedad",
y="Peso de granos por planta (g)"
)+
theme_bw()+
theme(
legend.position="none",
plot.title=element_text(face="bold",hjust=0.5),
axis.title=element_text(face="bold")
)
2.5 Distribución del rendimiento de los tratamientos por bloque
ggplot(datos,
aes(x = Tratamiento,
y = Peso,
fill = Tratamiento)) +
geom_boxplot(width = 0.65,
alpha = 0.8,
outlier.colour = "red") +
geom_jitter(width = 0.12,
size = 2,
alpha = 0.8) +
facet_wrap(~Humedad, nrow = 1) +
labs(
title = "Distribución del rendimiento por tratamiento en cada bloque de humedad",
subtitle = "Diseño en Bloques Completos al Azar",
x = "Tratamientos",
y = "Peso de granos por planta (g)"
) +
theme_bw(base_size = 13) +
theme(
legend.position = "none",
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
strip.background = element_rect(fill = "grey90"),
strip.text = element_text(face = "bold"),
panel.grid.minor = element_blank()
)
2.6 Rendimiento promedio por tratamiento
ggplot(descriptivos,
aes(x=Tratamiento,
y=Media,
fill=Tratamiento))+
geom_col()+
geom_text(aes(label=round(Media,2)),
vjust=-0.4,
size=4)+
labs(
title="Peso promedio por tratamiento",
x="Tratamientos",
y="Peso promedio (g)"
)+
theme_bw()+
theme(
legend.position="none",
plot.title=element_text(face="bold",hjust=0.5)
)
## Coeficiente de variación
CV_general <- (sd(datos$Peso)/mean(datos$Peso))*100
cat("Coeficiente de variación general:",round(CV_general,2),"%\n")## Coeficiente de variación general: 67.5 %if(CV_general<10){
cat("La variabilidad experimental es baja.\n")
}else if(CV_general<20){
cat("La variabilidad experimental es moderada.\n")
}else{
cat("La variabilidad experimental es alta.\n")
}## La variabilidad experimental es alta.3 Análisis de varianza (ANOVA)
3.1 Hipótesis
El análisis de varianza permite determinar si existen diferencias significativas entre los tratamientos evaluados.
3.1.1 Hipótesis para tratamientos
- H₀: Todos los tratamientos presentan el mismo rendimiento promedio.
- H₁: Al menos un tratamiento presenta un rendimiento promedio diferente.
Se rechaza H₀ cuando el valor p es menor que 0.05.
3.2 Ajuste del modelo
#---------------------------------------------------------
# Modelo DBCA
#---------------------------------------------------------
modelo <- aov(Peso ~ Tratamiento + Humedad,
data = datos)
summary(modelo)## Df Sum Sq Mean Sq F value Pr(>F)
## Tratamiento 2 99 49.4 0.179 0.8367
## Humedad 3 2858 952.8 3.458 0.0286 *
## Residuals 30 8267 275.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.3 Tabla ANOVA
anova_resultados <- summary(modelo)
knitr::kable(
anova_resultados[[1]],
digits = 4,
caption = "Análisis de varianza (DBCA)"
)| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Tratamiento | 2 | 98.8558 | 49.4279 | 0.1794 | 0.8367 |
| Humedad | 3 | 2858.4263 | 952.8088 | 3.4578 | 0.0286 |
| Residuals | 30 | 8266.5166 | 275.5506 | NA | NA |
3.4 Interpretación automática del efecto de tratamientos
## =====================================## EFECTO DE LOS TRATAMIENTOS## =====================================## Valor p = 0.8367if(p_trat < 0.05){
cat("Conclusión:\n")
cat("Se rechaza la hipótesis nula.\n")
cat("Existen diferencias significativas entre los tratamientos.\n")
}else{
cat("Conclusión:\n")
cat("No se rechaza la hipótesis nula.\n")
cat("No existen diferencias significativas entre los tratamientos.\n")
}## Conclusión:
## No se rechaza la hipótesis nula.
## No existen diferencias significativas entre los tratamientos.3.5 Evaluación de la eficiencia del bloqueo
El valor F de los bloques permite evaluar si el criterio de bloqueo contribuyó a explicar parte de la variabilidad experimental.
## =====================================## EFICIENCIA DEL BLOQUEO## =====================================## Valor F = 3.458if(F_bloque > 1){
cat("El valor F es mayor que 1.\n")
cat("El bloqueo fue eficiente, ya que explicó una proporción de la variabilidad mayor que el error experimental.\n")
}else{
cat("El valor F es menor o igual a 1.\n")
cat("El bloqueo no fue eficiente, ya que no logró reducir la variabilidad experimental respecto al error.\n")
}## El valor F es mayor que 1.
## El bloqueo fue eficiente, ya que explicó una proporción de la variabilidad mayor que el error experimental.4 Comparación de medias mediante la prueba de Tukey
Una vez detectadas diferencias significativas mediante el ANOVA, se realizó la prueba de comparación múltiple de Tukey con un nivel de significancia del 5 %, con el fin de identificar cuáles tratamientos difieren estadísticamente.
4.1 Hipótesis
- H₀: No existen diferencias entre las medias de los tratamientos comparados.
- H₁: Existen diferencias entre al menos dos medias.
4.2 Prueba de Tukey
#---------------------------------------------------------
# Prueba de Tukey
#---------------------------------------------------------
tukey <- HSD.test(modelo,
trt = "Tratamiento",
group = TRUE)
tukey## $statistics
## MSerror Df Mean CV MSD
## 275.5506 30 26.53028 62.56895 16.70666
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey Tratamiento 3 3.48642 0.05
##
## $means
## Peso std r se Min Max Q25 Q50 Q75
## 2,5% 28.27583 16.45773 12 4.791925 3.49 54.57 13.240 28.43 35.7550
## 5% 27.01167 19.78797 12 4.791925 3.08 72.41 12.170 26.99 34.0700
## Control 24.30333 18.67988 12 4.791925 6.34 77.07 13.915 19.65 31.3725
##
## $comparison
## NULL
##
## $groups
## Peso groups
## 2,5% 28.27583 a
## 5% 27.01167 a
## Control 24.30333 a
##
## attr(,"class")
## [1] "group"4.3 Grupos de comparación
knitr::kable(
tukey$groups,
digits = 3,
caption = "Agrupación de tratamientos según la prueba de Tukey"
)| Peso | groups | |
|---|---|---|
| 2,5% | 28.276 | a |
| 5% | 27.012 | a |
| Control | 24.303 | a |
4.4 Comparaciones múltiples y valores p
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Peso ~ Tratamiento + Humedad, data = datos)
##
## $Tratamiento
## diff lwr upr p adj
## 5%-2,5% -1.264167 -17.97083 15.44250 0.9810110
## Control-2,5% -3.972500 -20.67916 12.73416 0.8285165
## Control-5% -2.708333 -19.41500 13.99833 0.9159729knitr::kable(
as.data.frame(comparaciones$Tratamiento),
digits = 3,
caption = "Comparaciones múltiples de Tukey con valores p"
)| diff | lwr | upr | p adj | |
|---|---|---|---|---|
| 5%-2,5% | -1.264 | -17.971 | 15.442 | 0.981 |
| Control-2,5% | -3.972 | -20.679 | 12.734 | 0.829 |
| Control-5% | -2.708 | -19.415 | 13.998 | 0.916 |
4.5 Interpretación automática
## =====================================## PRUEBA DE TUKEY## =====================================## Los tratamientos que comparten la misma letra## no presentan diferencias estadísticas significativas.## Los tratamientos con letras diferentes## sí presentan diferencias significativas (α = 0.05).4.6 Gráfico de medias con letras de Tukey
medias <- tukey$groups
medias$Tratamiento <- rownames(medias)
ggplot(medias,
aes(x = reorder(Tratamiento, -Peso),
y = Peso,
fill = Tratamiento)) +
geom_col(width = 0.7) +
geom_text(aes(label = groups),
vjust = -0.6,
size = 6,
fontface = "bold") +
labs(
title = "Comparación de medias mediante la prueba de Tukey",
x = "Tratamientos",
y = "Peso promedio (g)"
) +
theme_bw() +
theme(
legend.position = "none",
plot.title = element_text(face = "bold", hjust = 0.5),
axis.title = element_text(face = "bold")
)
5 Verificación de supuestos
Antes de interpretar los resultados del ANOVA, es necesario comprobar que los residuos del modelo cumplen los supuestos de normalidad y homogeneidad de varianzas.
5.1 Normalidad de los residuos (Shapiro-Wilk)
5.1.1 Hipótesis
- H₀: Los residuos siguen una distribución normal.
- H₁: Los residuos no siguen una distribución normal.
Se rechaza H₀ cuando p < 0.05.
#---------------------------------------------------------
# Prueba de Shapiro-Wilk
#---------------------------------------------------------
shapiro <- shapiro.test(residuals(modelo))
shapiro##
## Shapiro-Wilk normality test
##
## data: residuals(modelo)
## W = 0.8957, p-value = 0.0026085.2 Curva de normalidad de los residuos
# Histograma de los residuos con densidad
hist(residuals(modelo),
probability = TRUE,
col = "lightblue",
border = "black",
main = "Distribución de los residuos",
xlab = "Residuos")
# Curva de densidad de los residuos
lines(density(residuals(modelo)),
col = "blue",
lwd = 2)
# Curva de la distribución normal teórica
curve(dnorm(x,
mean = mean(residuals(modelo)),
sd = sd(residuals(modelo))),
add = TRUE,
col = "red",
lwd = 2,
lty = 2)
legend("topright",
legend = c("Densidad observada",
"Distribución normal"),
col = c("blue","red"),
lwd = 2,
lty = c(1,2),
bty = "n")
## =====================================## PRUEBA DE SHAPIRO-WILK## =====================================## Valor p = 0.0026if(shapiro$p.value > 0.05){
cat("No se rechaza la hipótesis nula.\n")
cat("Los residuos presentan una distribución aproximadamente normal.\n")
}else{
cat("Se rechaza la hipótesis nula.\n")
cat("Los residuos no presentan una distribución normal.\n")
}## Se rechaza la hipótesis nula.
## Los residuos no presentan una distribución normal.5.3 Homogeneidad de varianzas (Bartlett)
5.3.1 Hipótesis
- H₀: Las varianzas son homogéneas.
- H₁: Al menos una varianza difiere.
##
## Bartlett test of homogeneity of variances
##
## data: Peso by Tratamiento
## Bartlett's K-squared = 0.36666, df = 2, p-value = 0.8325## =====================================## PRUEBA DE BARTLETT## =====================================## Valor p = 0.8325if(bartlett$p.value > 0.05){
cat("No se rechaza la hipótesis nula.\n")
cat("Las varianzas pueden considerarse homogéneas.\n")
}else{
cat("Se rechaza la hipótesis nula.\n")
cat("No existe homogeneidad de varianzas.\n")
}## No se rechaza la hipótesis nula.
## Las varianzas pueden considerarse homogéneas.5.4 Gráfico Q-Q
6 Conclusiones generales
## =====================================## CONCLUSIONES GENERALES## =====================================# Tratamientos
if(anova_resultados[[1]]$`Pr(>F)`[1] < 0.05){
cat("- Se encontraron diferencias significativas entre tratamientos.\n")
}else{
cat("- No se encontraron diferencias significativas entre tratamientos.\n")
}## - No se encontraron diferencias significativas entre tratamientos.# Normalidad
if(shapiro$p.value > 0.05){
cat("- El supuesto de normalidad fue satisfactorio.\n")
}else{
cat("- El supuesto de normalidad no se cumplió.\n")
}## - El supuesto de normalidad no se cumplió.# Homogeneidad
if(bartlett$p.value > 0.05){
cat("- Se verificó homogeneidad de varianzas.\n")
}else{
cat("- No se verificó homogeneidad de varianzas.\n")
}## - Se verificó homogeneidad de varianzas.# Bloques
if(F_bloque > 1){
cat("- El criterio de bloqueo fue eficiente para controlar parte de la variabilidad experimental.\n")
}else{
cat("- El criterio de bloqueo no fue eficiente.\n")
}## - El criterio de bloqueo fue eficiente para controlar parte de la variabilidad experimental.##
## Análisis estadístico finalizado correctamente.7 Comentario y conclusión final del experimento
Bajo las condiciones del presente experimento, los tratamientos evaluados no presentaron diferencias estadísticamente significativas respecto al control (p > 0.05). Por lo tanto, desde un punto de vista agronómico, no se recomienda la aplicación de ninguna de las dosis de caolín evaluadas, ya que no generan un incremento significativo en el rendimiento y su uso implicaría costos adicionales innecesarios.