# =============================================================================
# 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 : [Nombre del estudiante]
# Curso : Diseño y Análisis 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)
## Warning: package 'ggplot2' was built under R version 4.5.3
library(car)
## Warning: package 'car' was built under R version 4.5.3
## Cargando paquete requerido: carData
## Warning: package 'carData' was built under R version 4.5.3
library(multcompView)
## Warning: package 'multcompView' was built under R version 4.5.3
# =============================================================================
# 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
# =============================================================================
# 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)
# =============================================================================
# 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))
# =============================================================================
# 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.
# =============================================================================
# 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
# =============================================================================
# 10. CONCLUSIONES
# =============================================================================
# 1. Se cumplieron los tres supuestos del ANOVA (normalidad, homocedasticidad
# e independencia), por lo que el analisis de varianza es valido.
# 2. El bloqueo (variabilidad del lote) resulto significativo: justifica el
# uso de un DBCA en lugar de un diseno completamente al azar (DCA).
# 3. La dosis de fertilizante foliar Ca-B tuvo un efecto altamente
# significativo sobre el numero de frutos por planta de uchuva.
# 4. La prueba de Tukey permite identificar cuales dosis difieren entre si
# y, con ello, recomendar la dosis mas conveniente agronomica y
# economicamente (ver tabla de medias y agrupamiento resultante).
# =============================================================================
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