1 .- Introducción

El presente documento detalla el análisis estadístico de un experimento diseñado bajo un Diseño Completamente al Azar (DCA) con un número desigual de repeticiones (diseño desbalanceado). El objetivo principal es comparar los rendimientos (expresados en tn/ha) de cuatro variedades de aceitunas de mesa: Sevillana, Ascolana, Empeltre e Ilo.

El análisis de varianza (ANOVA) evalúa la siguiente hipótesis estadística:

  • Hipótesis Nula (\(H_0\)): \(\mu_{Sevillana} = \mu_{Ascolana} = \mu_{Empeltre} = \mu_{Ilo}\) (Los rendimientos medios son iguales poblacionalmente).
  • Hipótesis Alternativa (\(H_1\)): Al menos un \(\mu_i\) es diferente a los demás.

2 .- Diseño Experimental (Croquis)

Para garantizar la correcta ejecución del experimento en las 24 unidades experimentales (parcelas), se ha realizado una asignación completamente aleatoria de los tratamientos. El siguiente gráfico representa la distribución en el campo.

library(ggplot2)
library(dplyr)

# Definir los tratamientos
tratamientos <- c(rep("Sevillana", 6), rep("Ascolana", 7), rep("Empeltre", 5), rep("Ilo", 6))

# Aleatorizar
set.seed(123) 
croquis_aleatorio <- sample(tratamientos)

# Crear datos del croquis (4x6 = 24 parcelas)
datos_croquis <- data.frame(
  Fila = rep(1:4, each = 6),
  Columna = rep(1:6, times = 4),
  Variedad = croquis_aleatorio
)

# Gráfico del croquis
ggplot(datos_croquis, aes(x = Columna, y = Fila, fill = Variedad, label = Variedad)) +
  geom_tile(color = "white", linewidth = 2) +
  geom_text(color = "white", fontface = "bold", size = 4.5) +
  scale_fill_manual(values = c("Sevillana" = "#6E8B3D", "Ascolana" = "#4A708B", 
                               "Empeltre" = "#CD853F", "Ilo" = "#8B2323")) +
  scale_y_reverse() + 
  theme_void() +
  labs(title = "Croquis del Campo Experimental (DCA)",
       subtitle = "Asignación aleatoria de 24 parcelas",
       caption = "N: Sevillana (6), Ascolana (7), Empeltre (5), Ilo (6)") +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
    plot.subtitle = element_text(hjust = 0.5, size = 12, margin = margin(b=10)),
    plot.caption = element_text(hjust = 0.5, face = "italic"),
    legend.position = "none"
  )

3 .- Preparación de Datos y Herramientas

3.1 .- Carga de Librerías

Para garantizar un análisis robusto y una presentación gráfica y tabular de alta calidad académica, se emplean los siguientes paquetes:

library(car)         # Prueba de Levene 
library(lmtest)      # Prueba de Durbin-Watson 
library(knitr)       # Tablas
library(kableExtra)  # Formato avanzado de tablas
library(broom)       # Ordenar salidas estadísticas
library(gridExtra)   # Arreglo de gráficos

3.2 .- Ingreso de Datos Empíricos

Los datos experimentales provienen de las mediciones de campo. A continuación, se construye el conjunto de datos interno:

# Rendimientos (tn/ha)
rendimiento <- c(4.9, 8.5, 9.4, 9.8, 7.4, 10.5,      # Sevillana
                 8.2, 4.7, 5.6, 6.3, 6.4, 5.2, 4.8,  # Ascolana
                 5.6, 4.9, 5.0, 4.7, 4.3,            # Empeltre
                 10.6, 9.6, 8.7, 10.3, 9.9, 8.7)     # Ilo

# Factor "variedad"
variedad <- factor(c(rep("Sevillana", 6), rep("Ascolana", 7), 
                     rep("Empeltre", 5), rep("Ilo", 6)))

datos_aceitunas <- data.frame(variedad, rendimiento)

4 .- Análisis Exploratorio

Antes de proceder con la inferencia estadística, es fundamental caracterizar la muestra mediante estadística descriptiva para comprender el comportamiento inicial de las variedades.

resumen_datos <- datos_aceitunas %>%
  group_by(variedad) %>%
  summarise(
    `N` = n(),
    `Media (tn/ha)` = round(mean(rendimiento), 2),
    `Desviación Estándar` = round(sd(rendimiento), 2),
    `CV (%)` = round((sd(rendimiento)/mean(rendimiento))*100, 2)
  )

resumen_datos %>%
  kbl(caption = "Tabla 1. Resumen Estadístico del Rendimiento", align = "c") %>%
  kable_classic(full_width = FALSE, html_font = "Arial") %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2")
Tabla 1. Resumen Estadístico del Rendimiento
variedad N Media (tn/ha) Desviación Estándar CV (%)
Ascolana 7 5.89 1.22 20.73
Empeltre 5 4.90 0.47 9.68
Ilo 6 9.63 0.80 8.30
Sevillana 6 8.42 2.03 24.13

5 .- Inferencia Estadística (ANOVA)

Se ajusta un modelo de regresión lineal para realizar el Análisis de Varianza. Se asume un nivel de significancia de \(\alpha = 0.05\).

modelo_dca <- aov(rendimiento ~ variedad, data = datos_aceitunas)

tabla_anova <- tidy(modelo_dca) %>%
  mutate(
    p.value = ifelse(p.value < 0.001, "< 0.001", as.character(round(p.value, 4))),
    statistic = round(statistic, 2),
    sumsq = round(sumsq, 2),
    meansq = round(meansq, 2)
  )

colnames(tabla_anova) <- c("Fuente", "Grados Lib.", "Suma Cuad.", "Cuad. Medios", "Valor F", "Valor p")

tabla_anova %>%
  kbl(caption = "Tabla 2. Análisis de Varianza (ANOVA)", align = "c") %>%
  kable_classic(full_width = FALSE, html_font = "Arial") %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2") %>%
  footnote(general = "Si p < 0.05, se rechaza la hipótesis nula.")
Tabla 2. Análisis de Varianza (ANOVA)
Fuente Grados Lib. Suma Cuad. Cuad. Medios Valor F Valor p
variedad 3 82.89 27.63 16.42 < 0.001
Residuals 20 33.65 1.68 NA NA
Note:
Si p < 0.05, se rechaza la hipótesis nula.

6 .- Verificación de Supuestos del Modelo

La validez de los resultados paramétricos recae en el cumplimiento riguroso de tres supuestos metodológicos: normalidad, homocedasticidad e independencia.

6.1 Normalidad de los Residuos

Evaluada mediante la prueba de Shapiro-Wilk.

residuos <- residuals(modelo_dca)
shapiro_test <- shapiro.test(residuos)

data.frame(
  Estadístico_W = round(shapiro_test$statistic, 4),
  Valor_p = round(shapiro_test$p.value, 4)
) %>%
  kbl(caption = "Tabla 3. Prueba de Normalidad (Shapiro-Wilk)", align = "c") %>%
  kable_classic(full_width = FALSE, html_font = "Arial")
Tabla 3. Prueba de Normalidad (Shapiro-Wilk)
Estadístico_W Valor_p
W 0.9467 0.2299

6.2 .- Homocedasticidad (Igualdad de Varianzas)

Evaluada mediante la prueba robusta de Levene.

levene_test <- leveneTest(rendimiento ~ variedad, data = datos_aceitunas)

data.frame(
  Grados_Libertad = levene_test$Df[1],
  Estadístico_F = round(levene_test$`F value`[1], 4),
  Valor_p = round(levene_test$`Pr(>F)`[1], 4)
) %>%
  kbl(caption = "Tabla 4. Prueba de Homocedasticidad (Levene)", align = "c") %>%
  kable_classic(full_width = FALSE, html_font = "Arial")
Tabla 4. Prueba de Homocedasticidad (Levene)
Grados_Libertad Estadístico_F Valor_p
3 1.899 0.1623

6.3 .- Independencia de los Errores

Evaluada mediante la prueba de Durbin-Watson.

dw_test <- dwtest(lm(rendimiento ~ variedad, data = datos_aceitunas))

data.frame(
  Estadístico_DW = round(dw_test$statistic, 4),
  Valor_p = round(dw_test$p.value, 4)
) %>%
  kbl(caption = "Tabla 5. Prueba de Independencia (Durbin-Watson)", align = "c") %>%
  kable_classic(full_width = FALSE, html_font = "Arial")
Tabla 5. Prueba de Independencia (Durbin-Watson)
Estadístico_DW Valor_p
DW 1.6902 0.073

6.4 .- Diagnóstico Visual

Los análisis gráficos complementan las pruebas de hipótesis para corroborar la estructura de los residuos del modelo lineal ajustado.

datos_residuos <- data.frame(
  Ajustados = fitted(modelo_dca),
  Residuos = residuals(modelo_dca),
  Residuos_Std = rstandard(modelo_dca)
)

p1 <- ggplot(datos_residuos, aes(x = Ajustados, y = Residuos)) +
  geom_point(color = "#34495e", size = 2, alpha = 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "#e74c3c", linewidth = 1) +
  geom_smooth(method = "loess", se = FALSE, color = "#2980b9") +
  theme_classic() +
  labs(title = "Residuos vs. Ajustados", x = "Ajustados", y = "Residuos") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

p2 <- ggplot(datos_residuos, aes(sample = Residuos_Std)) +
  stat_qq(color = "#34495e", size = 2, alpha = 0.8) +
  stat_qq_line(color = "#e74c3c", linewidth = 1) +
  theme_classic() +
  labs(title = "Q-Q Normal", x = "Cuantiles Teóricos", y = "Residuos Estandarizados") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

grid.arrange(p1, p2, ncol = 2)

7 .- Análisis Post-Hoc (Tukey HSD)

Dado que el ANOVA resultó estadísticamente significativo, se ejecuta la Prueba de Rangos Múltiples de Tukey para identificar qué variedades difieren significativamente entre sí.

7.1 .- Diferencias Numéricas

tukey_test <- TukeyHSD(modelo_dca)

tukey_df <- as.data.frame(tukey_test$variedad) %>%
  mutate(
    `p adj` = ifelse(`p adj` < 0.001, "< 0.001", as.character(round(`p adj`, 4))),
    diff = round(diff, 3),
    lwr = round(lwr, 3),
    upr = round(upr, 3)
  )

tukey_df %>%
  kbl(caption = "Tabla 6. Comparaciones Múltiples (Tukey HSD)", align = "c") %>%
  kable_classic(full_width = FALSE, html_font = "Arial") %>%
  row_spec(which(as.numeric(tukey_test$variedad[, "p adj"]) < 0.05), bold = TRUE, color = "#e74c3c") %>%
  footnote(general = "Filas en rojo indican diferencias significativas (p < 0.05).")
Tabla 6. Comparaciones Múltiples (Tukey HSD)
diff lwr upr p adj
Empeltre-Ascolana -0.986 -3.112 1.140 0.5748
Ilo-Ascolana 3.748 1.728 5.767 < 0.001
Sevillana-Ascolana 2.531 0.511 4.551 0.011
Ilo-Empeltre 4.733 2.535 6.932 < 0.001
Sevillana-Empeltre 3.517 1.318 5.715 0.0012
Sevillana-Ilo -1.217 -3.313 0.879 0.388
Note:
Filas en rojo indican diferencias significativas (p < 0.05).

7.2 .- Visualización de Intervalos de Confianza

tukey_df_plot <- as.data.frame(tukey_test$variedad)
tukey_df_plot$Pares <- rownames(tukey_df_plot)
tukey_df_plot$Significativo <- ifelse(tukey_df_plot$`p adj` < 0.05, "Diferencia Significativa", "No Significativo")

ggplot(tukey_df_plot, aes(x = diff, y = reorder(Pares, diff), color = Significativo)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", linewidth = 0.8) +
  geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0.2, linewidth = 1) +
  geom_point(size = 3) +
  scale_color_manual(values = c("Diferencia Significativa" = "#e74c3c", "No Significativo" = "#7f8c8d")) +
  theme_classic(base_size = 12) +
  labs(
    title = "Intervalos de Confianza (95%) - Tukey HSD",
    x = "Diferencia en Medias (tn/ha)",
    y = "Pares Comparados"
  ) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "bottom",
    legend.title = element_blank()
  )

8 .- Conclusión Visual (Formato APA)

El siguiente gráfico de cajas sumariza la tendencia central y variabilidad del estudio, consolidando los hallazgos en un formato estandarizado para publicación.

ggplot(datos_aceitunas, aes(x = reorder(variedad, rendimiento, FUN = median), y = rendimiento)) +
  geom_boxplot(fill = "#fdfdfd", color = "black", outlier.shape = 19, outlier.size = 2, width = 0.5) +
  stat_summary(fun = mean, geom = "point", shape = 4, size = 3, color = "#c0392b", stroke = 1.2) + 
  labs(
    title = "Distribución del Rendimiento por Variedad",
    subtitle = "Las cruces rojas (x) señalan la media aritmética.",
    x = "Variedad de Aceituna",
    y = "Rendimiento (tn/ha)"
  ) +
  scale_y_continuous(limits = c(4, 11), breaks = seq(4, 12, by = 1)) +
  theme_classic(base_size = 13) +
  theme(
    text = element_text(family = "sans", color = "black"),
    plot.title = element_text(face = "bold", hjust = 0.5, size = 14),
    plot.subtitle = element_text(hjust = 0.5, size = 11, color = "grey40"),
    axis.title.x = element_text(margin = margin(t = 12), face = "bold"),
    axis.title.y = element_text(margin = margin(r = 12), face = "bold"),
    axis.text = element_text(color = "black", size = 12),
    axis.line = element_line(color = "black", linewidth = 0.5),
    axis.ticks = element_line(color = "black", linewidth = 0.5)
  )