Caso de Estudio: Nutrición y Producción de Leche Grado C en Panamá
Evolución de un DCA a un DBCA para Control de Varianza Biológica
Ing. Andrés Chang M.
2026-10-27
Carga de paquetes necesarios para todo el análisis
1.- Introducción y Planteamiento del Problema
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 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. Sin embargo, en las fincas Grado C panameñas, la respuesta productiva no solo depende de la dieta, sino de factores biológicos intrínsecos de hatos sumamente heterogéneos. Las vacas presentan variaciones extremas en su estado fisiológico, siendo el número de partos (etapa productiva) un factor crítico que define su potencial de síntesis láctea.
1.1.- Pregunta de Investigación (Enfoque Metodológico)
El presente caso de estudio tiene un objetivo dual. Primero, modelar estadísticamente la respuesta productiva de dietas basadas en gramíneas mejoradas y asociaciones silvopastoriles, utilizando datos simulados que reflejen los parámetros reportados por la literatura agropecuaria. Segundo, y más importante, demostrar metodológicamente cómo un análisis estadístico que ignora la varianza biológica (Diseño Completamente al Azar - DCA) fracasa frente a un modelo que controla el ruido ambiental aislando la heterogeneidad de los animales (Diseño de Bloques Completos al Azar - DBCA).
Para ello, diseñaremos un experimento virtual en bloques, estructurado en parcelas individuales de pastoreo con 20 vacas. Los tratamientos simulados son:
- T1: Hyparrhenia rufa (Control - Faragua).
- T2: Brachiaria brizantha cv. Toledo.
- T3: Panicum maximum cv. Mombasa.
- T4: Brachiaria spp. + Leguminosa arbustiva (Asociación silvopastoril).
2.- Diseño Experimental, Croquis de Campo y Simulación de Datos
Para que el ensayo refleje la realidad biológica de la finca, la investigación se establece físicamente bajo un Diseño de Bloques Completos al Azar (DBCA). Reconociendo que el hato es heterogéneo, agrupamos a las 20 vacas en 5 bloques biológicos estrictos según su número de partos (de 1 a 5 partos). Dentro de cada bloque, se sortean de manera aleatoria los cuatro tratamientos forrajeros.
2.1.- Generación del Croquis de Aleatorización
Utilizamos las funciones biométricas de la librería
agricolae para garantizar que la aleatorización cumpla con
las restricciones del DBCA. El siguiente croquis determina la identidad
exacta de los animales y las parcelas a recolectar en campo.
set.seed(123) # Semilla para reproducibilidad estricta
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 el diseño formal en bloques mediante agricolae
diseno_dbca <- design.rcbd(trt = niveles_trt, r = 5, seed = 123)
# Formateamos el libro de campo generado por agricolae usando setNames para mayor robustez
croquis_dbca <- diseno_dbca$book %>%
setNames(c("Parcela", "Bloque_Num", "Tratamiento_Forraje")) %>%
mutate(
Bloque_Etiqueta = factor(niveles_bloque[as.numeric(Bloque_Num)], levels = niveles_bloque),
Nombre_Plot = paste("Bloque", Bloque_Num, "\n(", Bloque_Etiqueta, ")")
) %>%
group_by(Bloque_Num) %>%
mutate(Posicion = row_number()) %>%
ungroup()
# Visualización del Croquis
ggplot(croquis_dbca, aes(x = factor(Posicion), y = forcats::fct_rev(Nombre_Plot), fill = Tratamiento_Forraje, label = Tratamiento_Forraje)) +
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 generada por software: 1 Tratamiento por Bloque",
x = "Posición Aleatoria de la Parcela en el Bloque",
y = "Identidad del Bloque Biológico (Vacas por Partos)"
) +
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 1. Mapa de Campo - Diseño de Bloques Completos al Azar (DBCA)
2.2.- Simulación Vinculada de la Matriz de Producción Lechera
A partir del croquis generado, procedemos a simular matemáticamente los resultados tras 60 días de pastoreo. Es fundamental destacar que la matriz de producción se construye directamente sobre el croquis aleatorizado, asegurando que los datos no estén desvinculados del diseño.
Integrando la literatura (ej. Aguilar-Pérez et al., 2019), parametrizamos los efectos teóricos para que los sistemas silvopastoriles y forrajes mejorados aporten mayor volumen lácteo. Así mismo, introduciremos el sesgo natural (varianza) provocado por el número de partos de la vaca.
# Efectos teóricos de la dieta forrajera (L/vaca/día) basados en literatura
efecto_trt <- c("T1_Faragua" = 4.0, "T2_Toledo" = 4.8, "T3_Mombasa" = 5.2, "T4_Asoc" = 6.0)
# Variación biológica introducida 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)
# Vinculamos estrictamente el diseño de campo de agricolae con los rendimientos
datos_exp <- croquis_dbca %>%
mutate(
Vaca_ID = paste0("V-", Parcela),
# Producción = Efecto Tratamiento + Efecto Bloque + Error Aleatorio
Produccion_L_dia = round(efecto_trt[Tratamiento_Forraje] + efecto_bloque[Bloque_Etiqueta] + rnorm(20, mean = 0, sd = 0.4), 2)
) %>%
select(Parcela, Vaca_ID, Partos = Bloque_Etiqueta, Tratamiento_Forraje, Produccion_L_dia)
# Matriz recolectada de campo lista para el análisis
datos_exp %>%
kbl(caption = "Tabla 1. Registro de Producción Lechera Simulada (L/día) vinculado al Croquis Aleatorio", align = "c") %>%
kable_classic(full_width = F, html_font = "Times New Roman") %>%
row_spec(0, bold = TRUE, background = "#f2f2f2")| Parcela | Vaca_ID | Partos | Tratamiento_Forraje | Produccion_L_dia |
|---|---|---|---|---|
| 101 | V-101 | 1_Parto | T3_Mombasa | 2.71 |
| 102 | V-102 | 1_Parto | T1_Faragua | 1.00 |
| 103 | V-103 | 1_Parto | T4_Asoc | 3.68 |
| 104 | V-104 | 1_Parto | T2_Toledo | 2.50 |
| 201 | V-201 | 2_Partos | T1_Faragua | 3.72 |
| 202 | V-202 | 2_Partos | T2_Toledo | 4.88 |
| 203 | V-203 | 2_Partos | T4_Asoc | 5.65 |
| 204 | V-204 | 2_Partos | T3_Mombasa | 4.61 |
| 301 | V-301 | 3_Partos | T4_Asoc | 7.76 |
| 302 | V-302 | 3_Partos | T2_Toledo | 5.96 |
| 303 | V-303 | 3_Partos | T3_Mombasa | 6.88 |
| 304 | V-304 | 3_Partos | T1_Faragua | 5.20 |
| 401 | V-401 | 4_Partos | T2_Toledo | 7.17 |
| 402 | V-402 | 4_Partos | T1_Faragua | 6.26 |
| 403 | V-403 | 4_Partos | T4_Asoc | 8.29 |
| 404 | V-404 | 4_Partos | T3_Mombasa | 7.86 |
| 501 | V-501 | 5_Partos | T1_Faragua | 2.83 |
| 502 | V-502 | 5_Partos | T4_Asoc | 5.10 |
| 503 | V-503 | 5_Partos | T2_Toledo | 4.24 |
| 504 | V-504 | 5_Partos | T3_Mombasa | 3.88 |
3.- Fase de Análisis 1: El Error de Modelar con DCA
Aunque el experimento se montó físicamente en bloques, un analista estadístico inexperto podría cometer el error metodológico de procesar la matriz ignorando el factor bloque. Al aplicar un modelo lineal simple de Diseño Completamente al Azar (DCA), se asume erróneamente que las 20 vacas son fisiológicamente idénticas. Todo el “ruido” biológico generado por la diferencia de partos se absorbe en el error residual del modelo matemático.
3.1.- Análisis de Varianza Paramétrico (DCA)
\(H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4\)
\(H_1: Al menos un Forraje es diferente\)
# Análisis de Varianza ignorando la variable biológica
modelo_dca <- aov(Produccion_L_dia ~ Tratamiento_Forraje, data = datos_exp)
tidy(modelo_dca) %>%
kbl(caption = "Tabla 2. ANOVA Ingenuo - Asumiendo erróneamente un 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 | 13.370 | 4.457 | 1.149 | 0.36 |
| Residuals | 16 | 62.057 | 3.879 | NA | NA |
Bajo este modelo matemático (\(p > 0.05\)), el investigador aceptaría que no existen evidencias para rechazar la hipótesis nula, concluyendo erróneamente que los forrajes mejorados no rinden más que el pasto Faragua. La enorme Suma de Cuadrados Residual delata que el modelo está ciego ante una fuente de variación crítica.
3.2.- Validación Gráfica y Formal de Supuestos (DCA)
Para evidenciar la deficiencia estructural de modelar esto como un DCA, procedemos a validar los supuestos de los residuos.
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
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))
shapiro_dca <- shapiro.test(residuals(modelo_dca))
levene_dca <- leveneTest(Produccion_L_dia ~ Tratamiento_Forraje, data = datos_exp)
data.frame(
Supuesto = c("Normalidad (Shapiro-Wilk)", "Homocedasticidad (Levene)"),
Estadístico = c(shapiro_dca$statistic, levene_dca$`F value`[1]),
Valor_p = c(shapiro_dca$p.value, levene_dca$`Pr(>F)`[1])
) %>%
kbl(caption = "Tabla 3. Pruebas de Supuestos para el modelo erróneo 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.9270 | 0.1350 |
| Homocedasticidad (Levene) | 0.0831 | 0.9682 |
La enorme dispersión vertical en el gráfico de residuos evidencia el problema: la varianza extrínseca no controlada (partos) destruye la sensibilidad de la prueba, a pesar de que formalmente pasen los supuestos estadísticos.
3.3.- Intento de Rescate: Kruskal-Wallis
El analista inexperto podría pensar que el problema de no hallar significancia se debe a una violación de normalidad (aunque el test formal la apruebe por poco) e intenta un método no paramétrico para muestras independientes.
kruskal_dca <- kruskal.test(Produccion_L_dia ~ Tratamiento_Forraje, data = datos_exp)
data.frame(
Prueba = "Kruskal-Wallis (Rangos)",
Chi_Cuadrado = kruskal_dca$statistic,
Valor_p = kruskal_dca$p.value
) %>%
kbl(caption = "Tabla 4. Alternativa No Paramétrica asumiendo independencia total", digits = 4, align = "c") %>%
kable_classic(full_width = F, html_font = "Times New Roman")| Prueba | Chi_Cuadrado | Valor_p | |
|---|---|---|---|
| Kruskal-Wallis chi-squared | Kruskal-Wallis (Rangos) | 2.4286 | 0.4883 |
El resultado sigue siendo no significativo. La matemática confirma que el problema no era distribucional, sino de especificación del modelo: se eligió un marco teórico biológicamente ciego para este tipo de hato ganadero.
4.- Fase de Análisis 2: Modelo Correcto (DBCA)
Para evaluar correctamente los mismos datos obtenidos en campo,
estructuramos el ANOVA respetando el diseño lógico original. Al integrar
la variable Partos como un bloque en la fórmula, el modelo
extrae el “ruido” que genera la diferencia de edad y lactancia,
permitiendo comparar los pastos de forma aislada y justa.
4.1.- Análisis de Varianza Paramétrico (DBCA)
# Análisis de Varianza integrando la covariable biológica (Bloque)
modelo_dbca <- aov(Produccion_L_dia ~ Tratamiento_Forraje + Partos, data = datos_exp)
tidy(modelo_dbca) %>%
kbl(caption = "Tabla 5. ANOVA Riguroso - Modelo 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 | 13.3697 | 4.4566 | 57.7369 | 0 |
| Partos | 4 | 61.1310 | 15.2828 | 197.9952 | 0 |
| Residuals | 12 | 0.9263 | 0.0772 | NA | NA |
El cambio metodológico es drástico. Al extraer de la Suma de Cuadrados Residual el inmenso peso aportado por los partos (\(p < 0.001\)), el Cuadrado Medio del Error cae drásticamente. Inmediatamente, la sensibilidad de la prueba se dispara, revelando que el Factor Forraje sí es altamente significativo (\(p < 0.001\)).
4.2.- Validación de Supuestos y Robustez del Modelo DBCA
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
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))
shapiro_dbca <- shapiro.test(residuals(modelo_dbca))
dw_dbca <- dwtest(modelo_dbca)
data.frame(
Supuesto = c("Normalidad (Shapiro-Wilk)", "Independencia (Durbin-Watson)"),
Estadístico = c(shapiro_dbca$statistic, dw_dbca$statistic),
Valor_p = c(shapiro_dbca$p.value, dw_dbca$p.value)
) %>%
kbl(caption = "Tabla 6. Verificación Formal de Supuestos - Residuos del 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.9368 | 0.2081 |
| DW | Independencia (Durbin-Watson) | 2.7471 | 0.8324 |
El gráfico de residuos muestra ahora un comportamiento compacto e impecable. Al cumplirse la normalidad y la independencia (\(p > 0.05\)), el modelo lineal por bloques es matemáticamente exacto.
4.3.- Prueba No Paramétrica Equivalente (Friedman)
Como respaldo metodológico y para darle robustez a la investigación, se aplica la prueba de Friedman. Esta alternativa no paramétrica está diseñada estrictamente para diseños en bloques y rangos, ofreciendo seguridad frente a cualquier posible anomalía en muestras biológicas pequeñas (Conover, 1999).
friedman_res <- friedman.test(Produccion_L_dia ~ Tratamiento_Forraje | Partos, data = datos_exp)
data.frame(
Prueba = "Friedman (Rangos en Bloques)",
Chi_Cuadrado = friedman_res$statistic,
Valor_p = friedman_res$p.value
) %>%
kbl(caption = "Tabla 7. Confirmación No Paramétrica del diseño DBCA", digits = 4, align = "c") %>%
kable_classic(full_width = F, html_font = "Times New Roman")| Prueba | Chi_Cuadrado | Valor_p | |
|---|---|---|---|
| Friedman chi-squared | Friedman (Rangos en Bloques) | 13.56 | 0.0036 |
El resultado (\(p < 0.05\)) ratifica que, controlando el efecto fisiológico de la vaca, las dietas no rinden lo mismo.
4.4.- Prueba Post-Hoc: Comparación de Tukey HSD
Habiendo validado rigurosamente el modelo paramétrico DBCA, aplicamos la prueba de Tukey para estructurar la jerarquía productiva que arrojó la simulación.
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. Jerarquía Productiva (Tukey HSD) resultante del Modelo Correcto", 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.10 | a |
| T3_Mombasa | 5.19 | b |
| T2_Toledo | 4.95 | b |
| T1_Faragua | 3.80 | c |
Los grupos definidos por Tukey evidencian que el modelo matemático logró aislar a la perfección las diferencias que configuramos al simular los datos: la asociación silvopastoril se ubica en el pico productivo, diferenciándose estadísticamente del control Faragua.
5.- Discusión y Conclusiones Metodológicas
La ejecución de este caso de estudio tuvo como eje central demostrar, a través de la simulación de datos, el profundo impacto que tiene la elección del modelo estadístico sobre la investigación agropecuaria. Las conclusiones derivadas abordan estrictamente la eficacia metodológica:
- Consecuencias de la Mala Especificación del Modelo (El Fracaso del DCA): Al ignorar la variabilidad fisiológica intrínseca (el número de partos de las vacas), la enorme dispersión infló la varianza residual, impidiendo que el ANOVA y la prueba de Kruskal-Wallis detectaran los incrementos productivos teóricos. Este hallazgo ilustra por qué muchas investigaciones de campo pueden fallar en probar sus hipótesis si el diseño estadístico no responde a la realidad biológica (Reuter & Moffet, 2016).
- Sensibilidad Analítica del Diseño en Bloques (DBCA): Al estructurar el análisis respetando el principio de los bloques, el modelo DBCA logró purificar el término de error. Extraer matemáticamente la varianza aportada por la edad productiva de las vacas permitió que el factor de interés (los pastos) brillara con alta significancia estadística (\(p < 0.001\)). Esto establece un precedente metodológico crítico: en fincas de Leche Grado C en Panamá, donde los hatos carecen de homogeneidad genética y fisiológica (Camarena, 2011; González, 2014), el uso de bloques no es una opción, ya que el bloqueo puede mejorar la precisión cuando existe heterogenidad y para garantizar el poder estadístico de la prueba.
- Confirmación Estadística de Postulados Teóricos: Los datos evaluados fueron simulados siguiendo fielmente los preceptos de la literatura regional (Aguilar-Pérez et al., 2019) sobre la superioridad productiva de las asociaciones silvopastoriles. El éxito del modelo DBCA, ratificado por las comprobaciones de supuestos formales, la prueba de Friedman y jerarquizado por la prueba post-hoc de Tukey, no busca descubrir empíricamente que un pasto es mejor que otro, sino validar que el armazón estadístico propuesto es capaz de medir y aislar eficazmente dichas mejoras en entornos biológicamente complejos.
6.- Referencias Bibliográficas
- Aguilar-Pérez, C., Solorio-Sánchez, F., Ku-Vera, J., Magaña-Monforte, J., & Santos-Flores, J. (2019). Producción de leche y carne en sistemas silvopastoriles. Bioagrociencias, 12(1), 1-8.
- Araúz, E. E. (2011). Alteración energética de deficiencia lactacional por estrés calórico. Curso de producción bovina de leche. Facultad de Ciencias Agropecuarias, Universidad de Panamá.
- Camarena D., M. I. (2011). Situación reproductiva y lactacional en fincas lecheras Grado C en los distritos de Alanje, Bugaba y Boquerón y su relación con la producción de leche y la condición bioeconómica. [Tesis de Licenciatura, Universidad de Panamá].
- Conover, W. J. (1999). Practical Nonparametric Statistics (3rd ed.). John Wiley & Sons.
- González N., M. M. (2014). Evaluación bioeconómica mediante un modelo de simulación (LIFE-SIM) para la producción bovina de leche en la zona baja de la Provincia de Chiriquí. [Tesis de Licenciatura, Universidad de Panamá].
- Montezuma, A. (2011). Diagnóstico y caracterización zootécnica de las fincas lecheras grado C en el distrito de Alanje. [Tesis de Licenciatura, Universidad de Panamá].
- Reuter, R. R., & Moffet, C. A. (2016). Invited Review: Designing a grazing experiment that can reliably detect meaningful differences. The Professional Animal Scientist, 32(1), 19-30.