1 Introducción

El Diseño de Bloques Completos al Azar (DBCA) —en inglés Randomized Complete Block Design (RCBD)— es uno de los diseños experimentales más usados en las ciencias agropecuarias. Su lógica es simple pero poderosa: cuando el área experimental (o el conjunto de animales) no es homogéneo, agrupamos las unidades experimentales en bloques internamente homogéneos y aplicamos todos los tratamientos dentro de cada bloque. Así, la variabilidad debida a la fuente de heterogeneidad conocida (bloque) se separa del error experimental, aumentando la precisión de las comparaciones entre tratamientos.

Modelo estadístico:

\[Y_{ij} = \mu + \tau_i + \beta_j + \varepsilon_{ij}\]

donde:

  • \(Y_{ij}\): respuesta observada del tratamiento \(i\) en el bloque \(j\),
  • \(\mu\): media general,
  • \(\tau_i\): efecto del tratamiento \(i\) (con \(i = 1, \dots, t\)),
  • \(\beta_j\): efecto del bloque \(j\) (con \(j = 1, \dots, b\)),
  • \(\varepsilon_{ij}\): error aleatorio, \(\varepsilon_{ij} \sim N(0, \sigma^2)\) independientes.

¿Cuándo bloquear? Cuando existe un factor de perturbación (nuisance) que podemos identificar y controlar pero que no nos interesa estudiar. Ejemplos típicos:

  • Agrícola: gradiente de fertilidad, humedad o pendiente del terreno; franjas de riego; fecha de siembra.
  • Pecuario: peso inicial del animal, camada/lote, corral, raza o edad.

Ejemplo de este documento (agrícola). Evaluaremos el rendimiento (t/ha) de cuatro híbridos de maíz (tratamientos) sembrados en un lote con un gradiente de fertilidad que dividimos en cinco bloques perpendiculares al gradiente.

Adaptación pecuaria. El mismo código sirve, por ejemplo, para evaluar la ganancia de peso (kg) de cuatro dietas en bovinos, bloqueando por peso inicial (bloque 1 = animales livianos, …, bloque 5 = pesados). Solo cambian los nombres de las variables y la interpretación.

2 Paquetes requeridos

# Ejecute este bloque una sola vez si no tiene los paquetes instalados.
install.packages(c("tidyverse", "agricolae", "car", "knitr"))
library(tidyverse)   # manipulación de datos y gráficos (dplyr, ggplot2)
library(agricolae)   # diseños experimentales y prueba de Duncan
library(car)         # prueba de Levene para homocedasticidad
library(knitr)       # tablas con kable()

3 Datos

Para que el documento sea reproducible, simulamos un conjunto de datos realista. En su curso puede reemplazar este bloque por la lectura de un archivo propio, por ejemplo:

datos <- read.csv("mis_datos.csv")
# datos <- readxl::read_excel("mis_datos.xlsx")
set.seed(2025)  # garantiza que todos obtengan los mismos números

tratamientos <- c("H1_ICA-305", "H2_Pioneer", "H3_Regional", "H4_Hibrido-X")
bloques       <- paste0("B", 1:5)

# Efectos "verdaderos" para la simulación (t/ha)
efecto_trat <- c("H1_ICA-305" = 0.0, "H2_Pioneer" = 1.2,
                 "H3_Regional" = 0.5, "H4_Hibrido-X" = 2.0)

# Gradiente de fertilidad entre bloques
efecto_bloq <- c("B1" = -0.8, "B2" = -0.4, "B3" = 0.0, "B4" = 0.4, "B5" = 0.8)

# Estructura completa: cada tratamiento aparece una vez por bloque
datos <- expand.grid(tratamiento = tratamientos,
                     bloque = bloques,
                     stringsAsFactors = FALSE)

media_general <- 7  # rendimiento base (t/ha)
datos$rendimiento <- media_general +
  efecto_trat[datos$tratamiento] +
  efecto_bloq[datos$bloque] +
  rnorm(nrow(datos), mean = 0, sd = 0.4)   # error experimental

datos$rendimiento <- round(datos$rendimiento, 2)

# Es CLAVE que tratamiento y bloque sean factores
datos$tratamiento <- factor(datos$tratamiento)
datos$bloque      <- factor(datos$bloque)

str(datos)
## 'data.frame':    20 obs. of  3 variables:
##  $ tratamiento: Factor w/ 4 levels "H1_ICA-305","H2_Pioneer",..: 1 2 3 4 1 2 3 4 1 2 ...
##  $ bloque     : Factor w/ 5 levels "B1","B2","B3",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ rendimiento: num  6.45 7.41 7.01 8.71 6.75 7.73 7.26 8.57 6.86 8.48 ...
##  - attr(*, "out.attrs")=List of 2
##   ..$ dim     : Named int [1:2] 4 5
##   .. ..- attr(*, "names")= chr [1:2] "tratamiento" "bloque"
##   ..$ dimnames:List of 2
##   .. ..$ tratamiento: chr [1:4] "tratamiento=H1_ICA-305" "tratamiento=H2_Pioneer" "tratamiento=H3_Regional" "tratamiento=H4_Hibrido-X"
##   .. ..$ bloque     : chr [1:5] "bloque=B1" "bloque=B2" "bloque=B3" "bloque=B4" ...

Estructura del experimento en formato de tabla de doble entrada (tratamientos × bloques):

tabla_ancha <- datos %>%
  pivot_wider(names_from = bloque, values_from = rendimiento)

kable(tabla_ancha, caption = "Rendimiento (t/ha) por tratamiento y bloque")
Rendimiento (t/ha) por tratamiento y bloque
tratamiento B1 B2 B3 B4 B5
H1_ICA-305 6.45 6.75 6.86 7.23 7.78
H2_Pioneer 7.41 7.73 8.48 8.91 9.17
H3_Regional 7.01 7.26 7.34 8.33 7.77
H4_Hibrido-X 8.71 8.57 8.30 9.32 10.77

4 Análisis exploratorio de datos (EDA)

Antes de modelar, describimos y entendemos los datos.

4.1 Estadísticos descriptivos por tratamiento

resumen_trat <- datos %>%
  group_by(tratamiento) %>%
  summarise(
    n      = n(),
    media  = mean(rendimiento),
    de     = sd(rendimiento),
    cv     = 100 * sd(rendimiento) / mean(rendimiento),
    minimo = min(rendimiento),
    maximo = max(rendimiento),
    .groups = "drop"
  ) %>%
  arrange(desc(media))

kable(resumen_trat, digits = 2,
      caption = "Estadísticos descriptivos por tratamiento (t/ha)")
Estadísticos descriptivos por tratamiento (t/ha)
tratamiento n media de cv minimo maximo
H4_Hibrido-X 5 9.13 0.99 10.82 8.30 10.77
H2_Pioneer 5 8.34 0.75 9.03 7.41 9.17
H3_Regional 5 7.54 0.52 6.88 7.01 8.33
H1_ICA-305 5 7.01 0.51 7.29 6.45 7.78

4.2 Estadísticos descriptivos por bloque

Un efecto de bloque marcado justifica el uso del DBCA sobre un diseño completamente al azar.

resumen_bloque <- datos %>%
  group_by(bloque) %>%
  summarise(media = mean(rendimiento),
            de = sd(rendimiento),
            .groups = "drop")

kable(resumen_bloque, digits = 2,
      caption = "Rendimiento promedio por bloque (t/ha)")
Rendimiento promedio por bloque (t/ha)
bloque media de
B1 7.40 0.96
B2 7.58 0.77
B3 7.74 0.77
B4 8.45 0.91
B5 8.87 1.43

5 Visualización

5.1 Distribución del rendimiento por tratamiento

ggplot(datos, aes(x = reorder(tratamiento, rendimiento, FUN = median),
                  y = rendimiento, fill = tratamiento)) +
  geom_boxplot(alpha = 0.7, outlier.shape = 21) +
  geom_jitter(width = 0.12, size = 2, alpha = 0.6) +
  labs(title = "Rendimiento por tratamiento",
       x = "Tratamiento (híbrido)", y = "Rendimiento (t/ha)") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "none")

5.2 Efecto de bloque

ggplot(datos, aes(x = bloque, y = rendimiento, fill = bloque)) +
  geom_boxplot(alpha = 0.7) +
  labs(title = "Rendimiento por bloque (gradiente de fertilidad)",
       x = "Bloque", y = "Rendimiento (t/ha)") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "none")

5.3 Gráfico de interacción tratamiento × bloque

Líneas aproximadamente paralelas indican que el ordenamiento de los tratamientos es consistente entre bloques (lo esperado en un DBCA sin interacción).

ggplot(datos, aes(x = bloque, y = rendimiento,
                  color = tratamiento, group = tratamiento)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  labs(title = "Interacción tratamiento × bloque",
       x = "Bloque", y = "Rendimiento (t/ha)", color = "Tratamiento") +
  theme_minimal(base_size = 13)

6 Modelación

Ajustamos el modelo del DBCA con aov(). El orden importa a nivel de escritura, pero al ser un diseño balanceado la suma de cuadrados es la misma:

modelo <- aov(rendimiento ~ tratamiento + bloque, data = datos)
anova_tabla <- summary(modelo)
anova_tabla
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## tratamiento  3 12.916   4.305  25.666 1.66e-05 ***
## bloque       4  6.283   1.571   9.364  0.00113 ** 
## Residuals   12  2.013   0.168                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretación del ANOVA. Si el valor \(p\) asociado a tratamiento es menor que \(\alpha = 0.05\), existe al menos un par de tratamientos con medias diferentes. El término bloque confirma si el bloqueo fue efectivo (su significancia respalda la decisión de bloquear).

Coeficiente de variación (CV) del experimento —indicador de precisión muy usado en agronomía:

cme <- summary(modelo)[[1]]["Residuals", "Mean Sq"]   # cuadrado medio del error
cv_exp <- 100 * sqrt(cme) / mean(datos$rendimiento)
cat("CME (error) =", round(cme, 4), "\n")
## CME (error) = 0.1677
cat("CV del experimento =", round(cv_exp, 2), "%\n")
## CV del experimento = 5.11 %
# cv = (sd/media)*100

En cultivos, un CV por debajo de ~15–20 % suele considerarse aceptable, aunque depende del cultivo y la variable medida.

7 Validación de supuestos

El ANOVA es válido si se cumplen: normalidad de los residuos, homocedasticidad (varianzas iguales) e independencia.

datos$ajustados <- fitted(modelo)
datos$residuos  <- residuals(modelo)

7.1 Normalidad de los residuos

Prueba de Shapiro-Wilk (H0: los residuos siguen distribución normal):

shapiro.test(datos$residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  datos$residuos
## W = 0.97149, p-value = 0.786

Gráfico Q-Q (los puntos deben alinearse sobre la recta):

ggplot(datos, aes(sample = residuos)) +
  stat_qq(color = "steelblue", size = 2) +
  stat_qq_line(color = "red", linewidth = 1) +
  labs(title = "Gráfico Q-Q de los residuos",
       x = "Cuantiles teóricos", y = "Cuantiles muestrales") +
  theme_minimal(base_size = 13)

7.2 Homocedasticidad

Prueba de Levene (robusta a la no normalidad) y de Bartlett (H0: las varianzas son iguales entre tratamientos):

leveneTest(rendimiento ~ tratamiento, data = datos)
bartlett.test(rendimiento ~ tratamiento, data = datos)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  rendimiento by tratamiento
## Bartlett's K-squared = 2.2404, df = 3, p-value = 0.524

Residuos vs. valores ajustados (la nube debe verse sin patrón ni forma de embudo):

ggplot(datos, aes(x = ajustados, y = residuos)) +
  geom_point(color = "steelblue", size = 2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Residuos vs. valores ajustados",
       x = "Valores ajustados (t/ha)", y = "Residuos") +
  theme_minimal(base_size = 13)

7.3 Independencia

En campo, la independencia se asegura sobre todo con la aleatorización de la asignación de tratamientos dentro de cada bloque. Como diagnóstico gráfico, revisamos los residuos en el orden de registro:

plot(datos$residuos, type = "b", pch = 19, col = "steelblue",
     xlab = "Orden de la unidad experimental", ylab = "Residuo",
     main = "Residuos en el orden de registro")
abline(h = 0, lty = 2, col = "red")

Si algún supuesto falla: considere transformar la variable respuesta (p. ej. \(\log\), raíz cuadrada), usar pruebas no paramétricas (Friedman para DBCA) o modelos mixtos.

8 Prueba de comparación de medias: Duncan

Cuando el ANOVA detecta diferencias significativas, la prueba de Duncan (Duncan’s Multiple Range Test) identifica cuáles tratamientos difieren entre sí. Es ampliamente usada en investigación agropecuaria por su sensibilidad para separar medias.

duncan_res <- duncan.test(modelo, "tratamiento",
                          group = TRUE, console = TRUE)
## 
## Study: modelo ~ "tratamiento"
## 
## Duncan's new multiple range test
## for rendimiento 
## 
## Mean Square Error:  0.1677475 
## 
## tratamiento,  means
## 
##              rendimiento       std r        se  Min   Max  Q25  Q50  Q75
## H1_ICA-305         7.014 0.5111066 5 0.1831652 6.45  7.78 6.75 6.86 7.23
## H2_Pioneer         8.340 0.7533923 5 0.1831652 7.41  9.17 7.73 8.48 8.91
## H3_Regional        7.542 0.5187196 5 0.1831652 7.01  8.33 7.26 7.34 7.77
## H4_Hibrido-X       9.134 0.9879929 5 0.1831652 8.30 10.77 8.57 8.71 9.32
## 
## Alpha: 0.05 ; DF Error: 12 
## 
## Critical Range
##         2         3         4 
## 0.5643882 0.5907525 0.6067262 
## 
## Means with the same letter are not significantly different.
## 
##              rendimiento groups
## H4_Hibrido-X       9.134      a
## H2_Pioneer         8.340      b
## H3_Regional        7.542      c
## H1_ICA-305         7.014      c

Tabla de agrupamiento (tratamientos que comparten letra no difieren estadísticamente al 5 %):

kable(duncan_res$groups, caption = "Prueba de Duncan: agrupamiento de medias")
Prueba de Duncan: agrupamiento de medias
rendimiento groups
H4_Hibrido-X 9.134 a
H2_Pioneer 8.340 b
H3_Regional 7.542 c
H1_ICA-305 7.014 c

Visualización del agrupamiento:

grupos <- duncan_res$groups %>%
  rownames_to_column("tratamiento") %>%
  rename(media = rendimiento, grupo = groups)

ggplot(grupos, aes(x = reorder(tratamiento, media), y = media, fill = grupo)) +
  geom_col(alpha = 0.85, width = 0.6) +
  geom_text(aes(label = paste0(round(media, 2), " (", grupo, ")")),
            hjust = -0.1, size = 4) +
  coord_flip() +
  labs(title = "Prueba de Duncan (α = 0.05)",
       subtitle = "Tratamientos con la misma letra no difieren significativamente",
       x = "Tratamiento", y = "Rendimiento medio (t/ha)", fill = "Grupo") +
  expand_limits(y = max(grupos$media) * 1.15) +
  theme_minimal(base_size = 13)

9 Conclusiones

A partir del análisis:

  1. El ANOVA permite decidir si los tratamientos (híbridos) difieren en rendimiento y si el bloqueo fue efectivo al controlar el gradiente de fertilidad.
  2. La validación de supuestos (normalidad, homocedasticidad e independencia) respalda la validez de las conclusiones del modelo.
  3. La prueba de Duncan ordena los tratamientos y agrupa aquellos estadísticamente equivalentes, orientando la recomendación práctica: seleccionar el híbrido (o dieta, en el caso pecuario) del grupo superior que además convenga por costo, disponibilidad o adaptación local.

Cierre para el aula. Reemplace la variable respuesta y los factores por los de su propio ensayo (agrícola o pecuario) y todo el flujo —EDA → visualización → modelo → supuestos → Duncan— se ejecuta sin modificaciones.

10 Referencias

  • Montgomery, D. C. Design and Analysis of Experiments. Wiley.
  • de Mendiburu, F. agricolae: Statistical Procedures for Agricultural Research. Paquete de R.
  • Kuehl, R. O. Diseño de experimentos: principios estadísticos de diseño y análisis de investigación.