=============================================================================

EFECTO DE LA FERTILIZACIÓN FOLIAR (Ca-B) SOBRE LA PRODUCCIÓN DE FRUTOS

EN UCHUVA (Physalis peruviana L.)

Diseño experimental : Bloques Completos al Azar (DBCA) - un solo factor

Elaborado por : Andres Santiago Suarez,Juan Manuel Anchique, Mateo Santiago Bustos, Juan David Millan

Curso : Diseño de experimentos

=============================================================================

# -----------------------------------------------------------------------------
# 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)

=============================================================================

1. DESCRIPCIÓN DEL EXPERIMENTO

=============================================================================

Pregunta de investigación:

¿La dosis de un fertilizante foliar a base de calcio-boro (Ca-B) afecta

el número de frutos producidos por planta de uchuva?

FACTOR DE ESTUDIO (1 factor): Dosis de fertilizante foliar Ca-B

NIVELES DEL FACTOR / TRATAMIENTOS (4 niveles):

T0 = 0.0 mL/L (testigo, se asperja solo agua)

T1 = 2.5 mL/L

T2 = 5.0 mL/L

T3 = 7.5 mL/L

FACTOR DE BLOQUEO (1 bloque): “Bloque” -> controla la variabilidad de

humedad y pendiente del lote experimental (5 niveles: B1…B5)

DISEÑO EXPERIMENTAL: Bloques Completos al Azar (DBCA), arreglo unifactorial

(un solo factor evaluado en varios niveles + 1 fuente de bloqueo).

Dentro de cada bloque, los 4 tratamientos se asignan al azar a una

parcela (aleatorización restringida por bloque).

REPETICIONES: 5 (una por bloque; cada bloque contiene los 4

tratamientos una sola vez)

UNIDADES EXPERIMENTALES: 4 tratamientos x 5 bloques = 20 parcelas

(cada parcela = promedio de 10 plantas evaluadas)

VARIABLE RESPUESTA: Número de frutos cosechados por planta

MODELO ESTADÍSTICO:

Y_ij = mu + tau_i + beta_j + epsilon_ij

Y_ij : número de frutos en el tratamiento i, bloque j

mu : media general

tau_i : efecto del i-ésimo nivel de dosis (i = 1,…,4)

beta_j : efecto del j-ésimo bloque (j = 1,…,5)

epsilon_ij: error experimental ~ N(0, sigma^2) iid

HIPÓTESIS (factor tratamiento):

H0: tau_1 = tau_2 = tau_3 = tau_4 = 0 (la dosis no afecta el N. frutos)

H1: al menos un tau_i != 0

HIPÓTESIS (factor de bloqueo):

H0: beta_1 = … = beta_5 = 0 (no hay efecto de bloque)

H1: al menos un beta_j != 0

Nota: al tratarse de UN solo factor (no un factorial de 2 o más factores

cruzados), el modelo NO incluye término de interacción. Se analizan de

forma directa el efecto del bloque y el efecto del factor tratamiento.

=============================================================================

=============================================================================

2. CONSTRUCCIÓN Y ALEATORIZACIÓN DEL DISEÑO DE CAMPO

=============================================================================

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

=============================================================================

3. SIMULACIÓN DE LA VARIABLE RESPUESTA

=============================================================================

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

=============================================================================

4. ESTADÍSTICA DESCRIPTIVA

=============================================================================

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.

=============================================================================

5. GRÁFICA DESCRIPTIVA (variable respuesta vs. tratamientos)

=============================================================================

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.

=============================================================================

6. AJUSTE DEL MODELO (necesario para obtener los residuales del ANOVA)

=============================================================================

modelo <- aov(frutos ~ bloque + tratamiento, data = datos)
residuos <- residuals(modelo)
ajustados <- fitted(modelo)

=============================================================================

7. VALIDACIÓN DE SUPUESTOS DEL ANOVA

=============================================================================

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.

=============================================================================

8. ANÁLISIS DE VARIANZA (ANOVA)

=============================================================================

Como los tres supuestos se cumplen, se interpreta directamente el ANOVA.

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.

=============================================================================

9. COMPARACIÓN MÚLTIPLE DE MEDIAS (Tukey HSD)

=============================================================================

Se realiza porque el factor tratamiento resultó significativo (p < 0.05).

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.

=============================================================================

10. CONCLUSION

=============================================================================

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