# -----------------------------------------------------------------------------
# 0. PAQUETES
# -----------------------------------------------------------------------------
# ggplot2 -> gráfica descriptiva profesional
# car -> prueba de Levene (homogeneidad de varianzas)
# multcompView -> letras de significancia (a, b, c) para la prueba de Tukey
paquetes <- c("ggplot2", "car", "multcompView")
faltantes <- paquetes[!paquetes %in% rownames(installed.packages())]
if (length(faltantes) > 0) install.packages(faltantes)
library(ggplot2)
library(car)
library(multcompView)
set.seed(2026)
dosis_mLL <- c(T0 = 0, T1 = 2.5, T2 = 5.0, T3 = 7.5)
tratamientos <- factor(names(dosis_mLL), levels = names(dosis_mLL))
id_bloques <- paste0("B", 1:5)
# Aleatorización de los tratamientos DENTRO de cada bloque (característica
# central del DBCA: cada bloque es una "repetición completa" del experimento)
campo <- do.call(rbind, lapply(id_bloques, function(b) {
data.frame(bloque = b,
tratamiento = sample(tratamientos, length(tratamientos)))
}))
campo$parcela <- 1:nrow(campo)
campo$bloque <- factor(campo$bloque, levels = id_bloques)
cat("=== Croquis de campo aleatorizado (DBCA) ===\n")
## === Croquis de campo aleatorizado (DBCA) ===
print(campo)
## bloque tratamiento parcela
## 1 B1 T0 1
## 2 B1 T3 2
## 3 B1 T2 3
## 4 B1 T1 4
## 5 B2 T0 5
## 6 B2 T2 6
## 7 B2 T3 7
## 8 B2 T1 8
## 9 B3 T3 9
## 10 B3 T0 10
## 11 B3 T2 11
## 12 B3 T1 12
## 13 B4 T1 13
## 14 B4 T3 14
## 15 B4 T2 15
## 16 B4 T0 16
## 17 B5 T2 17
## 18 B5 T0 18
## 19 B5 T1 19
## 20 B5 T3 20
Se definen efectos “verdaderos” (conocidos por el experimentador, ya que los datos son simulados) y un error experimental con desviación estándar pequeña, de modo que los supuestos del ANOVA se cumplan razonablemente.
mu <- 35 # media general
efecto_tratamiento <- c(T0 = 0, T1 = 8, T2 = 14, T3 = 16) # frutos adicionales
efecto_bloque <- c(B1 = -3, B2 = -1.5, B3 = 0, B4 = 1.5, B5 = 3)
sigma <- 2 # sd pequeña del error
campo$frutos <- with(campo,
round(mu +
efecto_tratamiento[as.character(tratamiento)] +
efecto_bloque[as.character(bloque)] +
rnorm(nrow(campo), mean = 0, sd = sigma))
)
datos <- campo[, c("bloque", "tratamiento", "frutos")]
cat("\n=== Base de datos del experimento ===\n")
##
## === Base de datos del experimento ===
print(datos)
## bloque tratamiento frutos
## 1 B1 T0 31
## 2 B1 T3 45
## 3 B1 T2 44
## 4 B1 T1 38
## 5 B2 T0 34
## 6 B2 T2 46
## 7 B2 T3 45
## 8 B2 T1 41
## 9 B3 T3 51
## 10 B3 T0 35
## 11 B3 T2 51
## 12 B3 T1 43
## 13 B4 T1 41
## 14 B4 T3 51
## 15 B4 T2 49
## 16 B4 T0 37
## 17 B5 T2 54
## 18 B5 T0 39
## 19 B5 T1 43
## 20 B5 T3 52
resumen <- aggregate(frutos ~ tratamiento, data = datos,
FUN = function(x) c(n = length(x), media = mean(x),
de = sd(x), cv = sd(x)/mean(x)*100))
resumen <- do.call(data.frame, resumen)
names(resumen) <- c("tratamiento", "n", "media", "de", "cv_pct")
cat("\n=== Resumen descriptivo por tratamiento ===\n")
##
## === Resumen descriptivo por tratamiento ===
print(resumen, row.names = FALSE)
## tratamiento n media de cv_pct
## T0 5 35.2 3.033150 8.616904
## T1 5 41.2 2.049390 4.974248
## T2 5 48.8 3.962323 8.119513
## T3 5 48.8 3.492850 7.157479
Se observa una tendencia creciente: al aumentar la dosis del fertilizante también aumenta el número promedio de frutos. Sin embargo, los tratamientos T2 y T3 presentan exactamente la misma media, lo que sugiere que incrementar la dosis de 5 a 7.5 mL/L no mejora la producción.
etiquetas <- c(T0 = "T0\n(0 mL/L)", T1 = "T1\n(2.5 mL/L)",
T2 = "T2\n(5.0 mL/L)", T3 = "T3\n(7.5 mL/L)")
g_desc <- ggplot(datos, aes(x = tratamiento, y = frutos, fill = tratamiento)) +
geom_boxplot(width = 0.55, alpha = 0.75, outlier.shape = NA) +
geom_jitter(aes(shape = bloque), width = 0.08, size = 2.6, color = "gray20") +
scale_x_discrete(labels = etiquetas) +
scale_fill_brewer(palette = "YlGn", guide = "none") +
labs(title = "Producción de frutos de uchuva según dosis de fertilizante foliar Ca-B",
subtitle = "Diseño de Bloques Completos al Azar (DBCA) - 5 bloques",
x = "Tratamiento (dosis de fertilizante foliar)",
y = "Número de frutos por planta",
shape = "Bloque") +
theme_minimal(base_size = 13) +
theme(plot.title = element_text(face = "bold", size = 12.5),
panel.grid.minor = element_blank())
print(g_desc)
ggsave("grafica_descriptiva_frutos.png", g_desc, width = 7.5, height = 5.5, dpi = 300)
Se aprecia un desplazamiento de las cajas hacia valores más altos conforme aumenta la dosis del fertilizante. Además, la dispersión es relativamente baja y no se observan valores atípicos importantes, indicando que los datos son consistentes entre bloques.
modelo <- aov(frutos ~ bloque + tratamiento, data = datos)
residuos <- residuals(modelo)
ajustados <- fitted(modelo)
cat("\n=========================================================\n")
##
## =========================================================
cat(" VALIDACIÓN DE SUPUESTOS \n")
## VALIDACIÓN DE SUPUESTOS
cat("=========================================================\n")
## =========================================================
## 7.1 Normalidad de los residuales (Shapiro-Wilk) --------------------------
prueba_normalidad <- shapiro.test(residuos)
cat("\n--- Supuesto 1: Normalidad de los residuales (Shapiro-Wilk) ---\n")
##
## --- Supuesto 1: Normalidad de los residuales (Shapiro-Wilk) ---
print(prueba_normalidad)
##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.94637, p-value = 0.3153
cat(ifelse(prueba_normalidad$p.value > 0.05,
"-> p > 0.05: NO se rechaza H0. Los residuales SI siguen una distribucion normal.\n",
"-> p <= 0.05: se rechaza H0. Los residuales NO siguen una distribucion normal.\n"))
## -> p > 0.05: NO se rechaza H0. Los residuales SI siguen una distribucion normal.
## 7.2 Homogeneidad de varianzas (Levene, sobre los residuales) -------------
prueba_homogeneidad <- leveneTest(residuos ~ tratamiento, data = datos)
cat("\n--- Supuesto 2: Homocedasticidad (prueba de Levene) ---\n")
##
## --- Supuesto 2: Homocedasticidad (prueba de Levene) ---
print(prueba_homogeneidad)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.1905 0.9013
## 16
p_levene <- prueba_homogeneidad$`Pr(>F)`[1]
cat(ifelse(p_levene > 0.05,
"-> p > 0.05: NO se rechaza H0. Las varianzas SI son homogeneas (homocedasticidad).\n",
"-> p <= 0.05: se rechaza H0. Las varianzas NO son homogeneas.\n"))
## -> p > 0.05: NO se rechaza H0. Las varianzas SI son homogeneas (homocedasticidad).
## 7.3 Independencia de los errores ------------------------------------------
cat("\n--- Supuesto 3: Independencia de los errores ---\n")
##
## --- Supuesto 3: Independencia de los errores ---
cat("Se garantiza por diseno: dentro de cada bloque, los tratamientos se\n")
## Se garantiza por diseno: dentro de cada bloque, los tratamientos se
cat("asignaron aleatoriamente a las parcelas (ver seccion 2). No existe\n")
## asignaron aleatoriamente a las parcelas (ver seccion 2). No existe
cat("relacion sistematica entre unidades experimentales.\n")
## relacion sistematica entre unidades experimentales.
cat("Se revisa ademas que no haya un patron en los residuales vs. el orden:\n")
## Se revisa ademas que no haya un patron en los residuales vs. el orden:
print(round(residuos, 2))
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## -0.2 0.2 -0.8 0.8 0.8 -0.8 -1.8 1.8 0.7 -1.7 0.7 0.3 -1.2 1.2 -0.8 0.8
## 17 18 19 20
## 1.7 0.3 -1.7 -0.3
# Panel grafico de diagnostico (normalidad, homogeneidad y valores atipicos)
png("diagnosticos_modelo.png", width = 900, height = 850, res = 130)
par(mfrow = c(2, 2))
plot(modelo)
dev.off()
## png
## 2
par(mfrow = c(1, 1))
Normalidad:Para evaluar la normalidad de los residuales se aplicó la prueba estadística de Shapiro-Wilk. El análisis arrojó un valor \(p\) de 0.315, y dado que este resultado es superior al nivel de significancia del 0.05, no se rechaza la hipótesis nula. Esto nos permite confirmar que los residuales presentan una distribución aproximadamente normal, cumpliendo a cabalidad con el primer supuesto.
Homogeneidad de varianzas:Este segundo supuesto se verificó utilizando la prueba de Levene. El resultado obtenido fue un valor \(p\) muy alto (0.901), lo que indica claramente que \(p > 0.05\) y, por lo tanto, las varianzas son estadísticamente homogéneas. En términos prácticos, esto significa que la variabilidad de los datos se mantiene similar entre todos los tratamientos evaluados, una condición indispensable para que el ANOVA sea válido.
Independencia:Finalmente, el supuesto de independencia no requiere validarse mediante una prueba estadística en el software, ya que es una condición que se asegura metodológicamente desde el trabajo en campo. Gracias a que se realizó una correcta aleatorización al asignar los tratamientos dentro de cada bloque, podemos garantizar que no existe dependencia ni sesgo sistemático entre las observaciones del experimento.
cat("\n=========================================================\n")
##
## =========================================================
cat(" ANALISIS DE VARIANZA (ANOVA) \n")
## ANALISIS DE VARIANZA (ANOVA)
cat("=========================================================\n")
## =========================================================
print(summary(modelo))
## Df Sum Sq Mean Sq F value Pr(>F)
## bloque 4 142.0 35.50 18.36 4.72e-05 ***
## tratamiento 3 651.8 217.27 112.38 4.76e-09 ***
## Residuals 12 23.2 1.93
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova_tabla <- summary(modelo)[[1]]
p_bloque <- anova_tabla["bloque", "Pr(>F)"]
p_trat <- anova_tabla["tratamiento", "Pr(>F)"]
cat("\n--- Interpretacion ---\n")
##
## --- Interpretacion ---
cat("Nota: al ser un diseno de UN solo factor + bloque, no existe termino\n")
## Nota: al ser un diseno de UN solo factor + bloque, no existe termino
cat("de interaccion que evaluar (eso solo aplica a arreglos factoriales de\n")
## de interaccion que evaluar (eso solo aplica a arreglos factoriales de
cat("2 o mas factores cruzados). Se interpretan bloque y tratamiento.\n\n")
## 2 o mas factores cruzados). Se interpretan bloque y tratamiento.
cat(ifelse(p_bloque < 0.05,
sprintf("BLOQUE: p = %.2e < 0.05 -> el bloqueo SI fue efectivo; hubo\n variabilidad significativa entre bloques, lo que confirma que\n controlar esta fuente de variacion (usar DBCA en vez de un diseno\n completamente al azar) fue una decision acertada.\n", p_bloque),
sprintf("BLOQUE: p = %.3f >= 0.05 -> no hubo diferencias significativas\n entre bloques.\n", p_bloque)))
## BLOQUE: p = 4.72e-05 < 0.05 -> el bloqueo SI fue efectivo; hubo
## variabilidad significativa entre bloques, lo que confirma que
## controlar esta fuente de variacion (usar DBCA en vez de un diseno
## completamente al azar) fue una decision acertada.
cat(ifelse(p_trat < 0.05,
sprintf("\nTRATAMIENTO: p = %.2e < 0.05 -> la dosis de fertilizante foliar\n Ca-B SI tuvo un efecto significativo sobre el numero de frutos por\n planta. Se procede con una prueba de comparacion multiple de medias.\n", p_trat),
sprintf("\nTRATAMIENTO: p = %.3f >= 0.05 -> no hay evidencia de efecto del\n factor tratamiento sobre el numero de frutos.\n", p_trat)))
##
## TRATAMIENTO: p = 4.76e-09 < 0.05 -> la dosis de fertilizante foliar
## Ca-B SI tuvo un efecto significativo sobre el numero de frutos por
## planta. Se procede con una prueba de comparacion multiple de medias.
Interpretación del factor Bloque:El análisis arrojó un valor \(p\) menor a 0.05 (\(p < 0.05\)), confirmando que el efecto del bloqueo fue estadísticamente significativo. Esto demuestra que en el terreno existía una variabilidad heterogénea entre sus zonas —como pendientes o gradientes de fertilidad—, por lo que la decisión de implementar un Diseño de Bloques Completos al Azar (DBCA) fue totalmente acertada para aislar ese ruido y evitar que alterara los resultados.
Interpretación del factor Tratamiento y Conclusión:Por su parte, la evaluación del tratamiento también obtuvo un valor \(p\) menor a 0.05 (\(p < 0.05\)), permitiéndonos rechazar con seguridad la hipótesis nula. Como conclusión principal, podemos afirmar con sustento estadístico que la dosis aplicada de fertilizante foliar a base de Calcio y Boro (Ca-B) sí influye de forma significativa en la producción final y el número de frutos por planta en el cultivo de uchuva.
if (p_trat < 0.05) {
cat("\n=========================================================\n")
cat(" PRUEBA DE TUKEY (comparaciones de a pares - tratamiento) \n")
cat("=========================================================\n")
tukey_trat <- TukeyHSD(modelo, "tratamiento")
print(tukey_trat)
png("tukey_tratamiento.png", width = 800, height = 600, res = 130)
plot(tukey_trat, las = 1, cex.axis = 0.8)
dev.off()
# Letras de significancia (estilo reporte agronomico): tratamientos que
# comparten letra NO difieren significativamente entre si (Tukey, 5%).
# Se reordenan para que "a" corresponda a la media mas alta (convencion
# habitual en reportes agronomicos, ej. paquete agricolae).
letras_crudo <- multcompLetters(tukey_trat$tratamiento[, "p adj"])$Letters
orden_medias <- names(sort(tapply(datos$frutos, datos$tratamiento, mean),
decreasing = TRUE))
secuencia <- unique(letras_crudo[orden_medias])
letras <- setNames(letters[match(letras_crudo, secuencia)], names(letras_crudo))
tabla_final <- merge(resumen, data.frame(tratamiento = names(letras),
letra = letras),
by = "tratamiento")
tabla_final <- tabla_final[order(-tabla_final$media), ]
cat("\n--- Tabla final de medias con letras de significancia (Tukey, 5%) ---\n")
print(tabla_final, row.names = FALSE)
}
##
## =========================================================
## PRUEBA DE TUKEY (comparaciones de a pares - tratamiento)
## =========================================================
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = frutos ~ bloque + tratamiento, data = datos)
##
## $tratamiento
## diff lwr upr p adj
## T1-T0 6.0 3.389167 8.610833 9.42e-05
## T2-T0 13.6 10.989167 16.210833 0.00e+00
## T3-T0 13.6 10.989167 16.210833 0.00e+00
## T2-T1 7.6 4.989167 10.210833 8.80e-06
## T3-T1 7.6 4.989167 10.210833 8.80e-06
## T3-T2 0.0 -2.610833 2.610833 1.00e+00
##
## --- Tabla final de medias con letras de significancia (Tukey, 5%) ---
## tratamiento n media de cv_pct letra
## T2 5 48.8 3.962323 8.119513 a
## T3 5 48.8 3.492850 7.157479 a
## T1 5 41.2 2.049390 4.974248 b
## T0 5 35.2 3.033150 8.616904 c
Análisis de la Prueba de Tukey: La prueba de comparación múltiple de Tukey nos permite categorizar el comportamiento agronómico de los tratamientos. En primer lugar, observamos que T2 y T3 comparten el mismo grupo estadístico (grupo a), lo que demuestra que no existen diferencias significativas entre aplicar 5.0 o 7.5 mL/L de fertilizante en términos de rendimiento.
En segundo lugar, ambos tratamientos (T2 y T3) lograron una producción de frutos significativamente superior a la de T1 y al testigo T0. Por último, aunque T1 tuvo un desempeño menor que las dosis más altas, logró superar de forma estadísticamente significativa al grupo de control (T0), confirmando que cualquier dosis de fertilización foliar es más efectiva que no aplicar el producto.
En conclusión, rechazamos la hipótesis nula, confirmando que el fertilizante foliar Ca-B incide significativamente en el rendimiento de la uchuva. La aplicación del diseño DBCA fue acertada para controlar la variabilidad del suelo, y tras verificar los supuestos del ANOVA, la prueba de Tukey demostró que las dosis de 5.0 mL/L y 7.5 mL/L logran el mayor número de frutos, sin diferencias estadísticas entre sí. Por lo tanto, la dosis de 5.0 mL/L se consolida como la alternativa más eficiente e insuperable en costo-beneficio, aunque la dosis baja (2.5 mL/L) podría evaluarse según las metas de producción de cada agricultor
cat("\n\nScript finalizado. Archivos generados:\n")
##
##
## Script finalizado. Archivos generados:
cat(" - grafica_descriptiva_frutos.png\n")
## - grafica_descriptiva_frutos.png
cat(" - diagnosticos_modelo.png\n")
## - diagnosticos_modelo.png
if (p_trat < 0.05) cat(" - tukey_tratamiento.png\n")
## - tukey_tratamiento.png