En los sistemas de producción bovina de doble propósito en Panamá (Leche Grado C), como los observados en los distritos de Alanje, Bugaba y Boquerón en la provincia de Chiriquí, la base de la alimentación es el pastoreo extensivo. La introducción de pasturas mejoradas (ej. Brachiaria brizantha y Panicum maximum) busca aumentar los rendimientos frente a pastos naturalizados como Hyparrhenia rufa (Faragua). La literatura señala que la calidad bromatológica de los pastos en el trópico disminuye drásticamente durante la época seca, limitando la expresión del potencial productivo (Araúz, 2011; Montezuma, 2011).
Asimismo, Aguilar-Pérez et al. (2019) destacan que los sistemas silvopastoriles y la inclusión de leguminosas forrajeras pueden mejorar los rendimientos lecheros al aumentar el aporte proteico, logrando producciones superiores de leche. Sin embargo, en las fincas de Leche Grado C en Panamá, la respuesta productiva no solo depende de la dieta (forraje), sino de factores biológicos intrínsecos de hatos heterogéneos (cruces de Brahman x Pardo Suizo x Holstein). En particular, las vacas presentan variaciones extremas en su estado fisiológico, siendo el número de partos (lactancias) un factor crítico que define su potencial de producción y requerimientos nutricionales.
¿El tipo de oferta forrajera (incluyendo gramíneas mejoradas y asociaciones con leguminosas) influye significativamente en el rendimiento lechero (L/vaca/día) en vacas de doble propósito tras un periodo de evaluación de 60 días en pastoreo? Además, ¿cómo la inclusión de la etapa productiva (número de partos agrupados en bloques) mejora la capacidad matemática del modelo estadístico para detectar dichas diferencias nutricionales frente a un Diseño Completamente al Azar (DCA)?
Para responder a esto, diseñamos un experimento con 20 vacas alojadas en parcelas individuales de pastoreo rotacional (para evitar pseudorreplicación y competencia). Los tratamientos forrajeros fueron:
Para llevar a cabo el estudio, se seleccionaron 20 vacas con distintos números de partos (de 1 a 5 partos). Tras un periodo de pastoreo de 60 días bajo las dietas asignadas, se registró la producción de leche promedio diaria.
A continuación, simulamos matemáticamente esta realidad de campo. En la simulación introduciremos intencionalmente una fuerte varianza biológica originada por la diferencia en la edad de las vacas (una vaca de 3 o 4 partos en su pico productivo dará mucha más leche que una novilla de primer parto o una vaca muy vieja, independientemente de lo que coma).
set.seed(123) # Semilla para reproducibilidad de la investigación
niveles_trt <- c("T1_Faragua", "T2_Toledo", "T3_Mombasa", "T4_Asoc")
niveles_bloque <- c("1_Parto", "2_Partos", "3_Partos", "4_Partos", "5_Partos")
# Generamos 5 grupos biológicos (bloques) y 4 vacas por grupo
bloques <- factor(rep(niveles_bloque, each = 4))
# Asignamos los 4 tratamientos asegurando que haya 1 tratamiento por bloque
tratamientos <- factor(rep(niveles_trt, times = 5))
# Efecto real de la dieta forrajera (L/vaca/día) basado en literatura
efecto_trt <- c("T1_Faragua" = 4.0, "T2_Toledo" = 4.8, "T3_Mombasa" = 5.2, "T4_Asoc" = 6.0)
# Ruido biológico enorme por etapa productiva (L/vaca/día)
efecto_bloque <- c("1_Parto" = -2.5, "2_Partos" = -0.5, "3_Partos" = 1.5, "4_Partos" = 2.5, "5_Partos" = -1.0)
# Producción de leche resultante al día 60 (Efecto Dieta + Efecto Vaca + Error aleatorio)
produccion_leche <- efecto_trt[tratamientos] + efecto_bloque[bloques] + rnorm(20, mean = 0, sd = 0.4)
# Consolidamos la base de datos recolectada en campo
datos_exp <- data.frame(
Vaca_ID = paste0("V-", 1:20),
Partos = bloques,
Tratamiento_Forraje = tratamientos,
Produccion_L_dia = round(produccion_leche, 2)
)
# Mostramos los datos tal cual llegaron de campo tras los 60 días
datos_exp %>%
kbl(caption = "Tabla 1. Registro de Producción Lechera (L/día) Tras 60 Días de Pastoreo", align = "c") %>%
kable_classic(full_width = F, html_font = "Times New Roman") %>%
row_spec(0, bold = TRUE, background = "#f2f2f2")| Vaca_ID | Partos | Tratamiento_Forraje | Produccion_L_dia |
|---|---|---|---|
| V-1 | 1_Parto | T1_Faragua | 1.28 |
| V-2 | 1_Parto | T2_Toledo | 2.21 |
| V-3 | 1_Parto | T3_Mombasa | 3.32 |
| V-4 | 1_Parto | T4_Asoc | 3.53 |
| V-5 | 2_Partos | T1_Faragua | 3.55 |
| V-6 | 2_Partos | T2_Toledo | 4.99 |
| V-7 | 2_Partos | T3_Mombasa | 4.88 |
| V-8 | 2_Partos | T4_Asoc | 4.99 |
| V-9 | 3_Partos | T1_Faragua | 5.23 |
| V-10 | 3_Partos | T2_Toledo | 6.12 |
| V-11 | 3_Partos | T3_Mombasa | 7.19 |
| V-12 | 3_Partos | T4_Asoc | 7.64 |
| V-13 | 4_Partos | T1_Faragua | 6.66 |
| V-14 | 4_Partos | T2_Toledo | 7.34 |
| V-15 | 4_Partos | T3_Mombasa | 7.48 |
| V-16 | 4_Partos | T4_Asoc | 9.21 |
| V-17 | 5_Partos | T1_Faragua | 3.20 |
| V-18 | 5_Partos | T2_Toledo | 3.01 |
| V-19 | 5_Partos | T3_Mombasa | 4.48 |
| V-20 | 5_Partos | T4_Asoc | 4.81 |
En esta primera fase, un investigador inexperto asume erróneamente que las 20 vacas son un rebaño homogéneo (ignorando el número de partos). Asigna los tratamientos al azar en todo el terreno y analiza bajo un modelo lineal simple DCA.
# Diseño DCA con agricolae
diseno_dca <- design.crd(trt = niveles_trt, r = 5, seed = 42)
# Formatear datos para el gráfico del croquis
croquis_dca <- diseno_dca$book %>%
setNames(c("plots", "r", "Tratamiento")) %>%
mutate(
Fila = rep(1:4, times = 5),
Columna = rep(1:5, each = 4)
)
# Gráfico del Croquis DCA
ggplot(croquis_dca, aes(x = factor(Columna), y = factor(Fila), fill = Tratamiento, label = Tratamiento)) +
geom_tile(color = "black", width = 0.95, height = 0.95, size = 0.5) +
geom_text(color = "white", fontface = "bold", size = 3.5) +
scale_fill_manual(values = c("T1_Faragua" = "#e74c3c", "T2_Toledo" = "#3498db",
"T3_Mombasa" = "#f39c12", "T4_Asoc" = "#2ecc71")) +
scale_y_discrete(limits = rev(levels(factor(croquis_dca$Fila)))) +
theme_void() +
labs(
title = "Croquis DCA en Campo",
subtitle = "Aleatorización total: Ignora la etapa productiva (Partos) de la vaca."
) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 14, family = "Times New Roman"),
plot.subtitle = element_text(hjust = 0.5, size = 12, family = "Times New Roman", margin=margin(b=10)),
legend.position = "none"
)Figura 1. Mapa de Campo - Diseño Completamente al Azar (DCA)
Explicación del Diseño DCA: En un Diseño Completamente al Azar (DCA), las unidades experimentales (las 20 parcelas con sus respectivas vacas) se consideran totalmente homogéneas. La asignación de cada uno de los 4 tratamientos forrajeros se hace mediante un sorteo puro, sin restricciones, de modo que cualquier tratamiento tiene la misma probabilidad de caer en cualquier parcela. El problema en campo es que si, por azar, el Tratamiento 1 le toca a todas las vacas viejas y el Tratamiento 4 le toca a las vacas en su pico de producción, los resultados de la leche estarán completamente sesgados.
# Análisis de Varianza (ANOVA) - DCA
modelo_dca <- aov(Produccion_L_dia ~ Tratamiento_Forraje, data = datos_exp)
tidy(modelo_dca) %>%
kbl(caption = "Tabla 2. Tabla ANOVA - Diseño Completamente al Azar (DCA)", align = "c", digits = 3) %>%
kable_classic(full_width = F, html_font = "Times New Roman")| term | df | sumsq | meansq | statistic | p.value |
|---|---|---|---|---|---|
| Tratamiento_Forraje | 3 | 11.923 | 3.974 | 0.916 | 0.455 |
| Residuals | 16 | 69.422 | 4.339 | NA | NA |
Al observar la Tabla 2, el valor p es superior al umbral de 0.05. Bajo el enfoque de un DCA, no existen diferencias significativas entre las dietas forrajeras. La alta varianza de los residuos indica que hay mucho “ruido” oculto en el experimento provocado por ignorar a las vacas.
Procedemos a evaluar si el modelo paramétrico simple cumple con las exigencias teóricas (Normalidad, Homocedasticidad e Independencia).
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
# Gráficos Diagnósticos
plot(fitted(modelo_dca), residuals(modelo_dca), pch = 19, col = "#2c3e50",
main = "Residuos vs Ajustados (DCA)", xlab = "Valores Ajustados", ylab = "Residuos")
abline(h = 0, col = "red", lty = 2)
qqnorm(residuals(modelo_dca), pch = 19, col = "#2c3e50", main = "Q-Q Plot (DCA)")
qqline(residuals(modelo_dca), col = "red", lwd = 2)par(mfrow = c(1, 1))
# Pruebas Formales
shapiro_dca <- shapiro.test(residuals(modelo_dca))
levene_dca <- leveneTest(Produccion_L_dia ~ Tratamiento_Forraje, data = datos_exp)
dw_dca <- dwtest(modelo_dca)
data.frame(
Supuesto = c("Normalidad (Shapiro-Wilk)", "Homocedasticidad (Levene)", "Independencia (Durbin-Watson)"),
Estadístico = c(shapiro_dca$statistic, levene_dca$`F value`[1], dw_dca$statistic),
Valor_p = c(shapiro_dca$p.value, levene_dca$`Pr(>F)`[1], dw_dca$p.value)
) %>%
kbl(caption = "Tabla 3. Pruebas Formales de Supuestos para el DCA", digits = 4, align = "c") %>%
kable_classic(full_width = F, html_font = "Times New Roman")| Supuesto | Estadístico | Valor_p | |
|---|---|---|---|
| W | Normalidad (Shapiro-Wilk) | 0.9312 | 0.1628 |
| Homocedasticidad (Levene) | 0.0619 | 0.9791 | |
| DW | Independencia (Durbin-Watson) | 0.4523 | 0.0001 |
Interpretación: Aunque los supuestos de independencia y normalidad podrían mantenerse al límite (al ser los datos generados de una matriz con ruido), visualmente en los gráficos se observa una dispersión masiva inaceptable. Al haber ignorado el efecto biológico de los partos, la varianza residual está fuertemente inflada, demostrando que el DCA es inadecuado para entidades biológicas heterogéneas.
Dado que el DCA falló en detectar diferencias y su error es masivo, confirmamos si evaluando las medianas mediante rangos se detecta algún efecto, utilizando la prueba de Kruskal-Wallis.
kruskal_dca <- kruskal.test(Produccion_L_dia ~ Tratamiento_Forraje, data = datos_exp)
data.frame(
Prueba = "Kruskal-Wallis",
Chi_Cuadrado = kruskal_dca$statistic,
GL = kruskal_dca$parameter,
Valor_p = kruskal_dca$p.value
) %>%
kbl(caption = "Tabla 4. Prueba No Paramétrica Kruskal-Wallis (Enfoque DCA)", digits = 4, align = "c") %>%
kable_classic(full_width = F, html_font = "Times New Roman")| Prueba | Chi_Cuadrado | GL | Valor_p | |
|---|---|---|---|---|
| Kruskal-Wallis chi-squared | Kruskal-Wallis | 2.1187 | 3 | 0.5481 |
Incluso bajo un enfoque no paramétrico riguroso (\(p > 0.05\)), se concluye que ignorar la etapa fisiológica de las vacas imposibilita estadísticamente encontrar el valor real de los forrajes mejorados.
Para rescatar la investigación, identificamos que el factor causante de la enorme dispersión en los residuos es la Etapa Productiva (Número de Partos). Al estructurar el análisis integrando esta variable biológica como factor de “Bloque”, evolucionamos metodológicamente a un DBCA.
En el nuevo diseño en campo, el experimento está estructurado estrictamente para que cada uno de los 5 grupos de partos (Bloques) pruebe los 4 forrajes simultáneamente.
# Diseño DBCA con agricolae
diseno_dbca <- design.rcbd(trt = niveles_trt, r = 5, seed = 123)
# Preparar datos para graficar croquis DBCA
croquis_dbca <- diseno_dbca$book %>%
setNames(c("plots", "Bloque", "Tratamiento")) %>%
group_by(Bloque) %>%
mutate(Columna = row_number()) %>%
ungroup() %>%
mutate(Nombre_Bloque = paste("Bloque", Bloque, "\n(", niveles_bloque[as.numeric(Bloque)], ")"))
# Gráfico del Croquis DBCA
ggplot(croquis_dbca, aes(x = factor(Columna), y = forcats::fct_rev(Nombre_Bloque), fill = Tratamiento, label = Tratamiento)) +
geom_tile(color = "black", width = 0.95, height = 0.9, size = 0.5) +
geom_text(color = "white", fontface = "bold", size = 3.5) +
scale_fill_manual(values = c("T1_Faragua" = "#e74c3c", "T2_Toledo" = "#3498db",
"T3_Mombasa" = "#f39c12", "T4_Asoc" = "#2ecc71")) +
theme_minimal() +
labs(
title = "Croquis DBCA en Campo",
subtitle = "Aleatorización restringida: Cada Tratamiento aparece una vez por Bloque.",
x = "Posición Aleatoria dentro del Bloque",
y = "Identidad del Bloque (Etapa Productiva)"
) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 14, family = "Times New Roman"),
plot.subtitle = element_text(hjust = 0.5, size = 12, family = "Times New Roman", margin=margin(b=10)),
axis.title = element_text(face = "bold", family = "Times New Roman"),
axis.text.y = element_text(face = "bold", size = 10, family = "Times New Roman"),
axis.text.x = element_blank(),
panel.grid = element_blank(),
legend.position = "bottom"
)Figura 2. Mapa de Campo - Diseño de Bloques Completos al Azar (DBCA)
Explicación del Diseño DBCA: En el Diseño de Bloques Completos al Azar, el investigador reconoce que las vacas no son idénticas. Agrupa a las vacas en “Bloques” homogéneos (ej. las 4 novillas de 1 parto juntas, las 4 vacas de 3 partos juntas, etc.). Luego, restringe el sorteo: dentro de cada bloque se sortean los 4 tratamientos para asegurar que todos los pastos sean probados de forma justa por vacas de la misma edad. Esto permite, en la fórmula matemática, aislar y restar el efecto del número de partos, purificando la comparación entre los forrajes.
# Ajustamos nuevamente el modelo lineal (ANOVA) integrando los Partos como Bloque
modelo_dbca <- aov(Produccion_L_dia ~ Tratamiento_Forraje + Partos, data = datos_exp)
tidy(modelo_dbca) %>%
kbl(caption = "Tabla 5. Tabla ANOVA - Diseño de Bloques Completos al Azar (DBCA)", align = "c", digits = 4) %>%
kable_classic(full_width = F, html_font = "Times New Roman") %>%
row_spec(1:2, bold = T, color = "black", background = "#e8f4f8")| term | df | sumsq | meansq | statistic | p.value |
|---|---|---|---|---|---|
| Tratamiento_Forraje | 3 | 11.9233 | 3.9744 | 20.3441 | 1e-04 |
| Partos | 4 | 67.0778 | 16.7695 | 85.8384 | 0e+00 |
| Residuals | 12 | 2.3443 | 0.1954 | NA | NA |
El resultado es categórico. Al remover la varianza que causaban los bloques (Partos, \(p < 0.001\)), el error residual disminuye drásticamente. Ahora, el factor Tratamiento de Forraje es altamente significativo (\(p < 0.001\)). De igual forma, se evidencia que existe una diferencia innegable en la producción de leche según la edad de la vaca (efecto de Bloque altamente significativo).
Antes de aplicar pruebas post-hoc, validamos las exigencias del modelo lineal paramétrico para garantizar que las conclusiones no tengan margen de error estadístico.
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
# Gráficos Diagnósticos
plot(fitted(modelo_dbca), residuals(modelo_dbca), pch = 19, col = "#2980b9",
main = "Residuos vs Ajustados (DBCA)", xlab = "Valores Ajustados", ylab = "Residuos")
abline(h = 0, col = "red", lty = 2)
qqnorm(residuals(modelo_dbca), pch = 19, col = "#2980b9", main = "Q-Q Plot (DBCA)")
qqline(residuals(modelo_dbca), col = "red", lwd = 2)par(mfrow = c(1, 1))
# Pruebas Formales
shapiro_dbca <- shapiro.test(residuals(modelo_dbca))
levene_dbca <- leveneTest(Produccion_L_dia ~ Tratamiento_Forraje, data = datos_exp)
dw_dbca <- dwtest(modelo_dbca)
data.frame(
Supuesto = c("Normalidad (Shapiro-Wilk)", "Homocedasticidad (Levene)", "Independencia (Durbin-Watson)"),
Estadístico = c(shapiro_dbca$statistic, levene_dbca$`F value`[1], dw_dbca$statistic),
Valor_p = c(shapiro_dbca$p.value, levene_dbca$`Pr(>F)`[1], dw_dbca$p.value)
) %>%
kbl(caption = "Tabla 6. Verificación Formal de Supuestos - ANOVA DBCA", digits = 4, align = "c") %>%
kable_classic(full_width = F, html_font = "Times New Roman")| Supuesto | Estadístico | Valor_p | |
|---|---|---|---|
| W | Normalidad (Shapiro-Wilk) | 0.9674 | 0.6998 |
| Homocedasticidad (Levene) | 0.0619 | 0.9791 | |
| DW | Independencia (Durbin-Watson) | 2.2140 | 0.3543 |
Al cumplirse a cabalidad todos los supuestos paramétricos (\(p > 0.05\)) y evidenciarse en los gráficos una distribución residual limpia y aleatoria alrededor de cero, confirmamos que el modelo DBCA es exacto y matemáticamente válido.
Aunque el ANOVA paramétrico validó exitosamente todos sus supuestos, la literatura contemporánea sobre biometría agrícola sugiere que, para blindar las conclusiones académicas en muestras biológicas in vivo, se debe corroborar el hallazgo con la alternativa no paramétrica homóloga al DBCA: la Prueba de Friedman (Conover, 1999).
friedman_res <- friedman.test(Produccion_L_dia ~ Tratamiento_Forraje | Partos, data = datos_exp)
data.frame(
Prueba = "Friedman (Controlando Bloques)",
Chi_Cuadrado = friedman_res$statistic,
GL = friedman_res$parameter,
Valor_p = friedman_res$p.value
) %>%
kbl(caption = "Tabla 7. Prueba No Paramétrica de Friedman (Validación del DBCA)", digits = 4, align = "c") %>%
kable_classic(full_width = F, html_font = "Times New Roman")| Prueba | Chi_Cuadrado | GL | Valor_p | |
|---|---|---|---|---|
| Friedman chi-squared | Friedman (Controlando Bloques) | 11.9388 | 3 | 0.0076 |
Esta prueba arroja una clara significancia (\(p < 0.05\)), lo que corrobora sólidamente los resultados de nuestro ANOVA DBCA, dándole total rigor publicable a la investigación al afirmar, bajo dos paradigmas estadísticos diferentes, que los forrajes sí afectan la producción lechera.
Para determinar con precisión la jerarquía de los forrajes y establecer cuáles son superiores, aplicamos la prueba de Diferencia Significativa Honesta de Tukey (Tukey HSD).
tukey_res <- HSD.test(modelo_dbca, "Tratamiento_Forraje", group = TRUE)
tukey_letras <- data.frame(
Forraje = rownames(tukey_res$groups),
Rendimiento_Promedio = round(tukey_res$groups$Produccion_L_dia, 2),
Grupo_Estadistico = tukey_res$groups$groups
)
tukey_letras %>%
kbl(caption = "Tabla 8. Prueba de Tukey HSD: Agrupación de Medias para Forrajes", align = "c") %>%
kable_classic(full_width = F, html_font = "Times New Roman") %>%
row_spec(1, bold=T, color="white", background="#2ca25f")| Forraje | Rendimiento_Promedio | Grupo_Estadistico |
|---|---|---|
| T4_Asoc | 6.04 | a |
| T3_Mombasa | 5.47 | ab |
| T2_Toledo | 4.73 | bc |
| T1_Faragua | 3.98 | c |
Interpretación de Tukey: Los tratamientos que no comparten la misma letra son estadísticamente diferentes entre sí al nivel del 5%. La tabla revela el claro dominio productivo de la asociación silvopastoril con leguminosas (Grupo a), ubicándose en el primer lugar de rendimiento lechero, superando de manera contundente al control de Faragua (Grupo c).
(Nota técnica: El factor bloque -Partos- también resultó significativo en el ANOVA. Esto corrobora biológicamente que las vacas en su tercer y cuarto parto produjeron más leche que las novillas o vacas añosas, validando la imperiosa necesidad de haber bloqueado este factor desde el croquis de campo).
La reestructuración metodológica de este caso de estudio permitió llegar a conclusiones firmes, publicables y directamente aplicables a los agronegocios del occidente de Panamá: