Se evalúa el rendimiento (t/ha) de 2 variedades de arroz:
expuestas a 3 niveles de fertilización nitrogenada: 0, 100 y 200 kg N/ha (N0, N100, N200), en un diseño de bloques completos al azar con estructura factorial (4 bloques, 24 unidades experimentales).
Preguntas guía de la sesión: ¿influye el nitrógeno en el rendimiento? ¿se comportan igual las dos variedades? ¿existe interacción entre variedad y nitrógeno (es decir, la respuesta al N depende de la variedad)?
# Ajusta la ruta si el archivo está en otra carpeta, p.ej. "datos/Datos_arroz.csv"
archivo <- "Datos_arroz.csv"
arroz <- read_delim(
archivo,
delim = ";",
locale = locale(decimal_mark = "."), # los decimales usan punto
escape_double = FALSE,
trim_ws = TRUE # limpia espacios sobrantes (p.ej. "N0 ")
)
# En lugar de View() (que NO funciona al tejer el documento) mostramos las
# primeras filas como tabla:
head(arroz) %>% kable(caption = "Primeras filas del conjunto de datos")
| Tratamiento | Variedad | Nitrogeno | Bloque | Rendimiento |
|---|---|---|---|---|
| T1 | VA | N0 | 1 | 7.5 |
| T1 | VA | N0 | 2 | 6.0 |
| T1 | VA | N0 | 3 | 7.0 |
| T1 | VA | N0 | 4 | 8.5 |
| T2 | VA | N100 | 1 | 8.5 |
| T2 | VA | N100 | 2 | 6.8 |
arroz <- arroz %>%
mutate(
Tratamiento = as.factor(Tratamiento),
Variedad = as.factor(Variedad),
# Orden explícito de los niveles de N (evita que queden alfabéticos):
Nitrogeno = factor(str_trim(Nitrogeno), levels = c("N0", "N100", "N200")),
Bloque = as.factor(Bloque)
)
glimpse(arroz)
## Rows: 24
## Columns: 5
## $ Tratamiento <fct> T1, T1, T1, T1, T2, T2, T2, T2, T3, T3, T3, T3, T4, T4, T4…
## $ Variedad <fct> VA, VA, VA, VA, VA, VA, VA, VA, VA, VA, VA, VA, VB, VB, VB…
## $ Nitrogeno <fct> N0, N0, N0, N0, N100, N100, N100, N100, N200, N200, N200, …
## $ Bloque <fct> 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4…
## $ Rendimiento <dbl> 7.5, 6.0, 7.0, 8.5, 8.5, 6.8, 7.3, 8.4, 7.6, 5.9, 7.3, 8.2…
resumen <- arroz %>%
group_by(Variedad, Nitrogeno) %>%
summarise(n = n(), Media = mean(Rendimiento), DE = sd(Rendimiento),
.groups = "drop")
resumen %>%
kable(digits = 2,
caption = "Rendimiento medio por combinación Variedad × Nitrógeno")
| Variedad | Nitrogeno | n | Media | DE |
|---|---|---|---|---|
| VA | N0 | 4 | 7.25 | 1.04 |
| VA | N100 | 4 | 7.75 | 0.83 |
| VA | N200 | 4 | 7.25 | 0.97 |
| VB | N0 | 4 | 6.50 | 0.91 |
| VB | N100 | 4 | 7.50 | 1.02 |
| VB | N200 | 4 | 8.00 | 1.27 |
El objetivo es revisar visualmente el efecto de cada factor sobre la respuesta antes de modelar.
# Efecto del nitrógeno
ggplot(arroz, aes(Nitrogeno, Rendimiento, fill = Nitrogeno)) +
geom_boxplot(alpha = 0.6) +
geom_jitter(width = 0.15, alpha = 0.7) +
labs(x = "Nitrógeno", y = "Rendimiento (t/ha)",
title = "Efecto del nitrógeno sobre el rendimiento") +
theme_minimal() + theme(legend.position = "none")
# Efecto de la variedad
ggplot(arroz, aes(Variedad, Rendimiento, fill = Variedad)) +
geom_boxplot(alpha = 0.6) +
geom_jitter(width = 0.15, alpha = 0.7) +
labs(x = "Variedad", y = "Rendimiento (t/ha)",
title = "Efecto de la variedad sobre el rendimiento") +
theme_minimal() + theme(legend.position = "none")
# Efecto combinado (los 6 tratamientos)
ggplot(arroz, aes(interaction(Variedad, Nitrogeno), Rendimiento, fill = Variedad)) +
geom_boxplot(alpha = 0.6) +
geom_jitter(width = 0.15, alpha = 0.7) +
labs(x = "Variedad : Nitrógeno", y = "Rendimiento (t/ha)",
title = "Efecto de la combinación Variedad × Nitrógeno") +
theme_minimal()
medias <- arroz %>%
group_by(Variedad, Nitrogeno) %>%
summarise(Media = mean(Rendimiento), .groups = "drop")
# Nitrógeno en el eje X, una línea por variedad
ggplot(medias, aes(Nitrogeno, Media, group = Variedad, color = Variedad)) +
geom_line(linewidth = 1) + geom_point(size = 2) +
labs(x = "Nitrógeno", y = "Rendimiento medio (t/ha)",
title = "Interacción: respuesta al nitrógeno según variedad") +
theme_minimal()
# Variedad en el eje X, una línea por nivel de N
ggplot(medias, aes(Variedad, Media, group = Nitrogeno, color = Nitrogeno)) +
geom_line(linewidth = 1) + geom_point(size = 2) +
labs(x = "Variedad", y = "Rendimiento medio (t/ha)",
title = "Interacción: diferencia entre variedades según nivel de N") +
theme_minimal()
Lectura: líneas paralelas ⇒ no hay interacción; líneas que se cruzan o divergen ⇒ hay interacción (la respuesta a un factor depende del otro).
modelo1 <- aov(Rendimiento ~ Bloque + Nitrogeno * Variedad, data = arroz)
summary(modelo1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Bloque 3 17.578 5.859 82.014 1.57e-09 ***
## Nitrogeno 2 3.000 1.500 20.995 4.49e-05 ***
## Variedad 1 0.042 0.042 0.583 0.456904
## Nitrogeno:Variedad 2 2.333 1.167 16.330 0.000172 ***
## Residuals 15 1.072 0.071
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Coeficiente de variación del experimento (calidad del ensayo):
cv.model(modelo1)
## [1] 3.624284
Un ANOVA solo es válido si los residuales cumplen normalidad, homogeneidad de varianzas e independencia.
# 1) Normalidad de los residuales
shapiro.test(residuals(modelo1))
##
## Shapiro-Wilk normality test
##
## data: residuals(modelo1)
## W = 0.95354, p-value = 0.3226
# 2) Homogeneidad de varianzas entre tratamientos
bartlett.test(Rendimiento ~ Tratamiento, data = arroz)
##
## Bartlett test of homogeneity of variances
##
## data: Rendimiento by Tratamiento
## Bartlett's K-squared = 0.55862, df = 5, p-value = 0.9898
# 3) Gráficos diagnósticos
par(mfrow = c(2, 2)); plot(modelo1); par(mfrow = c(1, 1))
Si p > 0.05 en Shapiro-Wilk y en Bartlett, los supuestos se cumplen y el ANOVA es confiable.
Regla de decisión:
Nitrogeno:Variedad es
significativa → se comparan las 6
combinaciones (efectos simples), no los factores por
separado.duncan.test(modelo1, "Nitrogeno", console = TRUE)
##
## Study: modelo1 ~ "Nitrogeno"
##
## Duncan's new multiple range test
## for Rendimiento
##
## Mean Square Error: 0.07144444
##
## Nitrogeno, means
##
## Rendimiento std r se Min Max Q25 Q50 Q75
## N0 6.875 0.9910312 8 0.09450162 5.5 8.5 6.000 7.00 7.500
## N100 7.625 0.8730079 8 0.09450162 6.1 8.5 7.175 7.70 8.425
## N200 7.625 1.1221790 8 0.09450162 5.9 9.3 7.050 7.75 8.275
##
## Alpha: 0.05 ; DF Error: 15
##
## Critical Range
## 2 3
## 0.2848586 0.2986087
##
## Means with the same letter are not significantly different.
##
## Rendimiento groups
## N100 7.625 a
## N200 7.625 a
## N0 6.875 b
duncan.test(modelo1, "Variedad", console = TRUE)
##
## Study: modelo1 ~ "Variedad"
##
## Duncan's new multiple range test
## for Rendimiento
##
## Mean Square Error: 0.07144444
##
## Variedad, means
##
## Rendimiento std r se Min Max Q25 Q50 Q75
## VA 7.416667 0.8973024 12 0.07716025 5.9 8.5 6.95 7.4 8.25
## VB 7.333333 1.1726685 12 0.07716025 5.5 9.3 6.25 7.5 8.05
##
## Alpha: 0.05 ; DF Error: 15
##
## Critical Range
## 2
## 0.232586
##
## Means with the same letter are not significantly different.
##
## Rendimiento groups
## VA 7.416667 a
## VB 7.333333 a
duncan_compara <- duncan.test(modelo1, c("Nitrogeno", "Variedad"), alpha = 0.05)
duncan_compara$groups |> rownames_to_column("Tratamiento")
maxy <- ceiling(max(duncan_compara$groups[[1]]) * 1.15)
bar.group(duncan_compara$groups, ylim = c(0, maxy),
xlab = "Tratamientos (Nitrógeno : Variedad)",
ylab = "Rendimiento medio (t/ha)", col = "purple")