El diseño completamente al azar (DCA) es el diseño experimental más simple para comparar \(k\) tratamientos o grupos. Las unidades experimentales se asignan aleatoriamente a los tratamientos, de modo que cada unidad tiene la misma probabilidad de recibir cualquier tratamiento.
Modelo estadístico
Para un factor con \(k\) niveles (tratamientos) y \(n_i\) réplicas en el tratamiento \(i\):
\[ y_{ij} = \mu + \tau_i + \varepsilon_{ij}, \quad i = 1, \dots, k; \quad j = 1, \dots, n_i \]
donde:
\(y_{ij}\) es la observación \(j\) del tratamiento \(i\),
\(\mu\) es la media global,
\(\tau_i\) es el efecto del tratamiento \(i\) (con \(\sum_{i=1}^k n_i \tau_i = 0\) si se impone restricción),
\(\varepsilon_{ij} \stackrel{\text{i.i.d.}}{\sim} N(0, \sigma^2)\) son los errores independientes.
Hipótesis en el ANOVA de un factor
Se desea contrastar si todos los efectos de tratamiento son nulos:
\[ H_0: \mu_1 = \mu_2 = \dots = \mu_k \quad \text{vs.} \quad H_1: \text{al menos un } \mu_i \text{ es diferente} \]
En términos de efectos: \(H_0: \tau_1 = \tau_2 = \dots = \tau_k = 0\).
Descomposición de la suma de cuadrados
La variabilidad total se particiona en variabilidad entre grupos y dentro de grupos:
\[SC_{\text{total}} = SC_{\text{tratamientos}} + SC_{\text{error}}\]
donde:
\[SC_{\text{total}} = \sum_{i=1}^k \sum_{j=1}^{n_i} (y_{ij} - \bar{y}_{\cdot\cdot})^2 = \sum_{i=1}^{k} \sum_{j=1}^{n_i} y_{ij}^2 - \frac{y_{\cdot \cdot}^2}{N} \quad \text{con} \quad gl = N - 1\]
\[SC_{\text{trat}} = \sum_{i=1}^k n_i (\bar{y}_{i\cdot} - \bar{y}_{\cdot\cdot})^2 \quad \text{con} \quad gl = k - 1\]
\[SC_{\text{error}} = \sum_{i=1}^k \sum_{j=1}^{n_i} (y_{ij} - \bar{y}_{i\cdot})^2 \quad \text{con} \quad gl = N - k\]
con \(N = \sum_{i=1}^k n_i\).
Estadístico de prueba
El estadístico \(F\) se construye como el cociente de los cuadrados medios:
\[ F = \frac{CM_{\text{trat}}}{CM_{\text{error}}} = \frac{SC_{\text{trat}}/(k-1)}{SC_{\text{error}}/(N-k)} \]
Bajo \(H_0\), \(F \sim F_{k-1,\, N-k}\).
Tabla ANOVA
| Fuente | Suma de cuadrados (SC) | Grados de libertad (gl) | Cuadrado medio (CM) | \(F\) |
|---|---|---|---|---|
| Tratamientos | \(SC_{\text{trat}}\) | \(k-1\) | \(CM_{\text{trat}}\) | \(F\) |
| Error | \(SC_{\text{error}}\) | \(N-k\) | \(CM_{\text{error}}\) | |
| Total | \(SC_{\text{total}}\) | \(N-1\) |
Regla de decisión (nivel \(\alpha\))
\[ \text{Rechazar } H_0 \text{ si } F > F_{\alpha,\, k-1,\, N-k} \]
Valor \(p\)
\[ p\text{-valor} = P(F_{k-1,\, N-k} > F_{\text{obs}} \mid H_0) \]
Ejemplo en R
Los datos que se presentan son rendimientos en toneladas por hectárea de un pasto con tres niveles de fertilización nitrogenada. El diseño fue completamente aleatorizado, con cinco repeticiones por tratamiento. A continuación se presenta la tabla con los niveles de nitrogeno y los rendimientos obtenidos.
| 1 | 2 | 3 |
|---|---|---|
| 14.823 | 25.151 | 32.605 |
| 14.676 | 25.401 | 32.460 |
| 14.720 | 25.131 | 32.256 |
| 14.514 | 25.031 | 32.669 |
| 15.065 | 25.267 | 32.111 |
# Cargar librerías necesarias
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.5.3library(dplyr)## Warning: package 'dplyr' was built under R version 4.5.3##
## Adjuntando el paquete: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggpubr)
library(RColorBrewer)
library(ggthemes)## Warning: package 'ggthemes' was built under R version 4.5.3library(patchwork)
# Crear el data frame con los datos
nitrogeno <- factor(rep(c(1, 2, 3), each = 5))
rendimiento <- c(14.823, 14.676, 14.720, 14.514, 15.065, # Nivel 1
25.151, 25.401, 25.131, 25.031, 25.267, # Nivel 2
32.605, 32.460, 32.256, 32.669, 32.111) # Nivel 3
# Crear data frame
datos <- data.frame(Nitrogeno = nitrogeno, Rendimiento = rendimiento)Estadísticas descriptivas
resumen_graf <- datos %>%
group_by(Nitrogeno) %>%
summarise(
n = n(),
Media = mean(Rendimiento),
Mediana = median(Rendimiento),
DE = sd(Rendimiento),
Error_std = DE / sqrt(n()),
IC_inf = Media - qt(0.975, n() - 1) * Error_std,
IC_sup = Media + qt(0.975, n() - 1) * Error_std,
Min = min(Rendimiento),
Max = max(Rendimiento),
.groups = "drop"
)
print(resumen_graf)## # A tibble: 3 × 10
## Nitrogeno n Media Mediana DE Error_std IC_inf IC_sup Min Max
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 5 14.8 14.7 0.204 0.0911 14.5 15.0 14.5 15.1
## 2 2 5 25.2 25.2 0.142 0.0635 25.0 25.4 25.0 25.4
## 3 3 5 32.4 32.5 0.235 0.105 32.1 32.7 32.1 32.7# Media global
media_global <- mean(rendimiento)
cat("\nMedia global:", round(media_global, 4), "\n")##
## Media global: 24.1253Análisis visual
# Paleta de colores personalizada
colores <- c("#1B9E77", "#D95F02", "#7570B3")
colores_suaves <- c("#A6CEE3", "#B2DF8A", "#FB9A99")
colores_gradiente <- c("#E41A1C", "#377EB8", "#4DAF4A")
# Gráfico Boxplot con violín y puntos individuales
ggplot(datos, aes(x = Nitrogeno, y = Rendimiento, fill = Nitrogeno)) +
geom_violin(alpha = 0.5, trim = FALSE, show.legend = FALSE) +
geom_boxplot(width = 0.2, alpha = 0.8, show.legend = FALSE,
outlier.color = "red", outlier.size = 3, outlier.shape = 8) +
geom_jitter(width = 0.15, size = 2.5, alpha = 0.6, color = "black") +
stat_summary(fun = mean, geom = "point", shape = 18, size = 4,
color = "darkred", show.legend = FALSE) +
scale_fill_manual(values = colores) +
labs(
title = "Distribución del rendimiento por nivel de nitrógeno",
subtitle = "Boxplot con violín y puntos individuales | Media indicada en rojo",
x = "Nivel de nitrógeno",
y = "Rendimiento (ton/ha)",
caption = "Fuente: Datos experimentales"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 10, hjust = 0.5, color = "gray50"),
plot.caption = element_text(size = 8, hjust = 0, color = "gray50")
)
# Gráfico Medias con intervalos de confianza
ggplot(resumen_graf, aes(x = Nitrogeno, y = Media, color = Nitrogeno)) +
geom_point(size = 5, stroke = 1.5) +
geom_errorbar(aes(ymin = IC_inf, ymax = IC_sup), width = 0.2, size = 1) +
geom_line(aes(group = 1), size = 0.8, linetype = "dashed", color = "gray50") +
scale_color_manual(values = colores_gradiente) +
labs(
title = "Evolución del rendimiento medio por nivel de nitrógeno",
subtitle = "Barras = Intervalos de confianza al 95%",
x = "Nivel de nitrógeno",
y = "Rendimiento medio (ton/ha)"
) +
theme_minimal() +
theme(legend.position = "none")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ANOVA
# Modelo ANOVA
modelo <- aov(Rendimiento ~ Nitrogeno, data = datos)
# Obtener tabla ANOVA
anova_result <- summary(modelo)
# Imprimir tabla ANOVA
print(anova_result)## Df Sum Sq Mean Sq F value Pr(>F)
## Nitrogeno 2 788.3 394.2 10131 <2e-16 ***
## Residuals 12 0.5 0.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Verificación de los supuestos del modelo ANOVA de un factor
El modelo ANOVA de un factor se basa en tres supuestos fundamentales:
- Independencia de las observaciones
- Normalidad de los residuos
- Homocedasticidad (igualdad de varianzas entre grupos)
La violación de estos supuestos puede afectar la validez de las conclusiones del análisis.
Verificación de la independencia
La independencia depende principalmente del diseño experimental:
Asignación aleatoria de las unidades experimentales a los tratamientos
Ausencia de dependencia temporal o espacial entre observaciones
Cada observación proviene de una unidad experimental diferente
Verificación de la normalidad de los residuos
La normalidad de los residuos es uno de los supuestos fundamentales del modelo ANOVA. Los residuos (\(\varepsilon_{ij}\)) deben seguir una distribución normal con media cero y varianza constante.
Hipótesis sobre los residuos
\[ \varepsilon_{ij} \sim N(0, \sigma^2) \]
Hipótesis de la prueba de normalidad:
\[ H_0: \text{Los residuos siguen una distribución normal} \] \[ H_1: \text{Los residuos no siguen una distribución normal} \]
Métodos gráficos para evaluar normalidad
El gráfico \(Q-Q\) compara los cuantiles teóricos de una distribución normal con los cuantiles observados de los residuos. Si los puntos se alinean aproximadamente sobre la línea de referencia, los residuos son normales.
# Prueba de normalidad global
shapiro_test <- shapiro.test(residuals(modelo))
print(shapiro_test)##
## Shapiro-Wilk normality test
##
## data: residuals(modelo)
## W = 0.97219, p-value = 0.8891# Prueba de normalidad por grupo
for(i in 1:3) {
cat("\nNivel", i, ":\n")
print(shapiro.test(rendimiento[nitrogeno == i]))
}##
## Nivel 1 :
##
## Shapiro-Wilk normality test
##
## data: rendimiento[nitrogeno == i]
## W = 0.97108, p-value = 0.8821
##
##
## Nivel 2 :
##
## Shapiro-Wilk normality test
##
## data: rendimiento[nitrogeno == i]
## W = 0.96469, p-value = 0.8402
##
##
## Nivel 3 :
##
## Shapiro-Wilk normality test
##
## data: rendimiento[nitrogeno == i]
## W = 0.94058, p-value = 0.6701# Gráficos de normalidad
# Q-Q plot
qqnorm(residuals(modelo),
main = "Gráfico Q-Q de residuos",
pch = 16,
col = "blue")
qqline(residuals(modelo),
col = "red",
lwd = 2)
# Q-Q plot con ggplot
ggplot(data.frame(residuos = residuals(modelo)), aes(sample = residuos)) +
stat_qq(size = 2, color = "blue") +
stat_qq_line(color = "red", size = 1) +
labs(title = "Q-Q Plot de residuos",
x = "Cuantiles teóricos",
y = "Cuantiles muestrales") +
theme_minimal()
# Histograma de residuos
hist(residuals(modelo),
breaks = "FD",
col = "lightblue",
border = "white",
main = "Histograma de residuos",
xlab = "Residuos",
probability = TRUE)
# Superponer curva de densidad normal teórica
curve(dnorm(x,
mean = mean(residuals(modelo)),
sd = sd(residuals(modelo))),
add = TRUE,
col = "red",
lwd = 2)
# Superponer densidad empírica
lines(density(residuals(modelo)),
col = "darkgreen",
lwd = 2,
lty = 2)
# Añadir leyenda
legend("topright",
legend = c("Normal teórica", "Densidad empírica"),
col = c("red", "darkgreen"),
lwd = 2,
lty = c(1, 2))
Verificación de la homocedasticidad (igualdad de varianzas)
La homocedasticidad es el supuesto de que las varianzas son iguales en todos los grupos de tratamiento. En el modelo ANOVA, se asume que:
\[ \sigma_1^2 = \sigma_2^2 = \dots = \sigma_k^2 = \sigma^2 \]
La violación de este supuesto (heterocedasticidad) puede afectar la tasa de error tipo I y la potencia del test F.
Hipótesis sobre la igualdad de varianzas
\[ H_0: \sigma_1^2 = \sigma_2^2 = \dots = \sigma_k^2 \] \[ H_1: \text{Al menos una varianza es diferente} \]
Métodos gráficos para evaluar homocedasticidad
El gráfico de residuos vs. valores ajustados diagnóstico más importante para detectar heterocedasticidad. La dispersión de los residuos debe ser constante a lo largo de todos los valores ajustados.
# Prueba de Levene (recomendada)
library(car)## Cargando paquete requerido: carData##
## Adjuntando el paquete: 'car'## The following object is masked from 'package:dplyr':
##
## recodelevene_test <- leveneTest(Rendimiento ~ Nitrogeno, data = datos, center = median)
print(levene_test)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.5372 0.5978
## 12# Gráfico de residuos vs. valores ajustados
plot(fitted(modelo), residuals(modelo),
xlab = "Valores ajustados (medias predichas)",
ylab = "Residuos",
main = "Residuos vs. Valores ajustados",
pch = 16,
col = as.numeric(nitrogeno), # CORREGIDO: variedat -> nitrogeno
cex = 1.2)
abline(h = 0, lty = 2, col = "red", lwd = 2)
# Añadir línea de tendencia suavizada
lines(lowess(fitted(modelo), residuals(modelo)),
col = "blue", lwd = 2)
# Boxplot de la variable respuesta por tratamiento
boxplot(rendimiento ~ nitrogeno, # CORREGIDO: variedad -> nitrogeno
col = c("lightblue", "lightgreen", "lightcoral"),
main = "Distribución de la variable respuesta por tratamiento",
xlab = "Tratamiento",
ylab = "Rendimiento",
notch = FALSE,
varwidth = FALSE)
# Boxplot de residuos por grupo
boxplot(residuals(modelo) ~ nitrogeno, # CORREGIDO: variedad -> nitrogeno
col = "lightyellow",
main = "Residuos por tratamiento",
xlab = "Tratamiento",
ylab = "Residuos",
border = "darkblue")
abline(h = 0, lty = 2, col = "red", lwd = 2)
# Boxplot con ancho proporcional a la raíz del tamaño muestral
boxplot(rendimiento ~ nitrogeno, # CORREGIDO: variedad -> nitrogeno
varwidth = TRUE,
col = "lightgreen",
main = "Boxplot con ancho proporcional a √n",
xlab = "Tratamiento",
ylab = "Rendimiento")
Comparaciones o pruebas de rango múltiples
Cuando el análisis de varianza (ANOVA) rechaza la hipótesis nula de igualdad de medias, únicamente sabemos que al menos un par de medias es diferente, pero no sabemos cuál(es) par(es) específicamente. Para identificar qué tratamientos difieren entre sí se utilizan las pruebas de comparaciones múltiples (o pruebas post hoc).
Principales pruebas de rango múltiple
| Prueba | Característica | Control de error |
|---|---|---|
| Tukey HSD | Compara todos los pares de medias | Error familiar (FWER) al nivel \(\alpha\) |
| Bonferroni | Ajusta \(\alpha\) dividiendo por número de comparaciones | FWER conservador |
| Scheffé | Permite contrastes complejos (combinaciones lineales) | FWER muy conservador |
| Fisher LSD | Solo se aplica si ANOVA es significativo | Error por comparación (no controla FWER) |
| Duncan | Procedimiento jerárquico por rangos | Error por comparación |
| Student–Newman–Keuls (SNK) | Similar a Tukey pero menos conservador | Error por comparación |
FWER (Family-Wise Error Rate): probabilidad de cometer al menos un error tipo I en el conjunto de todas las comparaciones.
Prueba de Tukey HSD (Honestly Significant Difference)
Es la prueba más utilizada cuando se comparan todos los pares posibles de medias.
Estadístico de prueba:
\[ \text{HSD} = q_{\alpha,\, k,\, N-k} \cdot \sqrt{\frac{CM_{\text{error}}}{2} \left( \frac{1}{n_i} + \frac{1}{n_j} \right)} \]
donde:
\(q_{\alpha,\, k,\, N-k}\) es el valor crítico de la distribución de rango estudentizado
\(k\) = número de tratamientos
\(N-k\) = grados de libertad del error
\(CM_{\text{error}}\) = cuadrado medio del error
\(n_i, n_j\) = tamaños de muestra de los grupos comparados
Regla de decisión:
\[ \text{Rechazar } H_0: \mu_i = \mu_j \quad \text{si} \quad |\bar{y}_{i\cdot} - \bar{y}_{j\cdot}| > \text{HSD} \]
Valor p ajustado (para cada comparación):
\[ p_{\text{ajustado}} = P\left(Q_{k,\, N-k} > \frac{|\bar{y}_{i\cdot} - \bar{y}_{j\cdot}|}{\sqrt{CM_{\text{error}} \left( \frac{1}{n_i} + \frac{1}{n_j} \right)/2}} \right) \]
tukey_result <- TukeyHSD(modelo, conf.level = 0.95)
cat("\n========== PRUEBA POST HOC (Tukey HSD) ==========\n")##
## ========== PRUEBA POST HOC (Tukey HSD) ==========print(tukey_result)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Rendimiento ~ Nitrogeno, data = datos)
##
## $Nitrogeno
## diff lwr upr p adj
## 2-1 10.4366 10.103773 10.769427 0
## 3-1 17.6606 17.327773 17.993427 0
## 3-2 7.2240 6.891173 7.556827 0# Gráfico de Tukey
tukey_df <- as.data.frame(tukey_result$Nitrogeno)
tukey_df$Comparacion <- rownames(tukey_df)
tukey_df$Significativo <- tukey_df$`p adj` < 0.05
ggplot(tukey_df, aes(x = Comparacion, y = diff, color = Significativo)) +
geom_point(size = 4) +
geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.2, size = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
scale_color_manual(values = c("TRUE" = "darkgreen", "FALSE" = "red"),
labels = c("TRUE" = "Significativo (p < 0.05)",
"FALSE" = "No significativo")) +
labs(
title = "Intervalos de confianza de Tukey (95%)",
subtitle = "Diferencias entre niveles de nitrógeno",
x = "Comparación",
y = "Diferencia de medias",
color = ""
) +
theme_minimal() +
coord_flip()
Corrección de Bonferroni
Divide el nivel de significancia global \(\alpha\) entre el número total de comparaciones \(m = \binom{k}{2}\).
Valor crítico:
\[ t_{\text{crit}} = t_{\alpha/(2m),\, N-k} \]
Intervalo de confianza simultáneo (nivel \(1-\alpha\)):
\[ (\bar{y}_{i\cdot} - \bar{y}_{j\cdot}) \pm t_{\alpha/(2m),\, N-k} \cdot \sqrt{CM_{\text{error}}\left(\frac{1}{n_i} + \frac{1}{n_j}\right)} \]
Relación con Bonferroni:
\[ p_{\text{ajustado}} = \min(m \cdot p_{\text{individual}}, \ 1) \]
# Comparaciones con Bonferroni
bonferroni_result <- pairwise.t.test(rendimiento, nitrogeno,
p.adjust.method = "bonferroni")
print(bonferroni_result)##
## Pairwise comparisons using t tests with pooled SD
##
## data: rendimiento and nitrogeno
##
## 1 2
## 2 < 2e-16 -
## 3 < 2e-16 1.4e-15
##
## P value adjustment method: bonferroniF_calculado <- anova_result[[1]]$`F value`[1]
p_valor <- anova_result[[1]]$`Pr(>F)`[1]
# Crear etiquetas de significancia según Tukey
letras <- c("a", "b", "c")
resumen_graf$Letra <- letras
p_final <- ggplot(resumen_graf, aes(x = Nitrogeno, y = Media, fill = Nitrogeno)) +
geom_col(width = 0.7, alpha = 0.85) +
geom_errorbar(aes(ymin = IC_inf, ymax = IC_sup), width = 0.2, size = 1, color = "black") +
geom_text(aes(label = Letra, y = IC_sup + 0.8), size = 7, fontface = "bold") +
scale_fill_manual(values = colores) +
labs(
title = "Rendimiento del pasto por nivel de fertilización nitrogenada",
subtitle = paste("ANOVA: F =", round(F_calculado, 2),
"| p =", ifelse(p_valor < 0.001, "< 0.001", round(p_valor, 4)),
"| Tukey: medias con letras diferentes difieren significativamente (α = 0.05)"),
x = "Nivel de nitrógeno",
y = "Rendimiento medio (ton/ha)",
caption = "Barras de error: IC 95%"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11, hjust = 0.5, color = "gray30"),
plot.caption = element_text(size = 9, hjust = 0),
legend.position = "none",
axis.title = element_text(face = "bold"),
panel.grid.major.x = element_blank()
)
print(p_final)